Showing error for fix shake in lammps

LAMMPS LAMMPS Beginners
1

Dear LAMMPS Community,

I am currently trying to perform Langevin dynamics simulation for rigid dimers using the version 29-Aug2024. I first equilibrate the dimers with very stiff spring (without the command fix shake or fix rigid) and use that equilibrated configuration as input for the production run. In the production run, I use the command fix shake, but I am getting an error. Is it because the dimer is linear molecule which throws error. However, I have seen examples in lammps for dimers where the dimers are rigid by fix shake. But for my input script, I am getting error for fix shake. Alternatively, if I use fix rigid with langevin, then it is not showing any error, but the time step is much smaller in that case. Is it possible to use fix shake with larger time step for the rigid dimers. Your help on this matter would be greatly appreciated. Below, I have attached the lammps file and the initial configuration.

input.data (3.5 MB)
Production.in (905 Bytes)

2

what error?

3

Thanks for your reply.

Here is the error file. It seems the error comes from the linear geometry of the dimer. So, I increased the spring constant from 30 to 100. But, still getting the same error.

error.lammps (6.1 KB)

4

It would have been much nicer, if you had just quoted the error (WARNING: Shake determinant < 0.0) instead of forcing me to first download and look at the file.

That cannot be. A dimer is by its very definition linear.

When using fix shake, the force constant for the bonded interaction is ignored. So that change has no effect.

It is not obvious to me what is triggering this issue. Perhaps @sjplimp has an idea and can comment. It would be extremely helpful for debugging, if you could reproduce this issue with a much smaller system, say with just some 10s or at the most a few 100 atoms and post that input deck here.

5

Thanks for your reply.

I did not use spring constant 100 for the production run where I am using the shake command. I am using spring constant 100 for the equilibration without using the shake command. I am sorry for that confusion.

For 100 atoms (or 50 rigid dimers), I don’t get any error. The simulation works well. In the simulation box, the dimers are not interacting and the potential energy is zero till time step 10000000. However, the error still persists only for 16000 rigid dimers .

Here, I am attaching the scripts which I have used for the equilibration and the production run.

Equilibration.in (1.3 KB)
input.data (10.4 KB)
mymol.txt (166 Bytes)
Production.in (905 Bytes)

6

Dear @akohlmey and @sjplimp

Can you please help for the error “WARNING: Shake determinant < 0.0” using fix shake for 16000 rigid dimers? For smaller systems like 100 dimers, there is no problem for running the simulations. Do you have any idea how to resolve that issue for larger systems like 16000 rigid dimers using fix shake?

7

In the meantime, you can try using bond_style zero instead of bond_style harmonic for your production run. Apart from the (very small) chance that that’s where the error is coming from (for some reason I can’t imagine), this would save you (a tiny little bit of) computational time.

8

I don’t think that the size itself has a direct impact on this situation.

The SHAKE determinant error is due to the system of constraint equations not converging, which often is due to too large forces (varying by more than 15 orders of magnitude) or pathological cases like having three atoms in a row where the force to move the center atom so that the angle is less than 180 degrees is infinite.

What may have an impact is that for a larger system with more atoms, this situation is more likely to happen, when you create atoms/molecules on random positions. So what you should try out is the “overlap” keyword of create_atoms (i would first set it to something like 0.3) which will avoid any close contacts in the initial geometry.

With random placement, there is a small but finite chance you could have by chance two molecules placed in a way that their bonds cross and since you are doing a 2d system, they cannot unoverlap in this situation.

9

I can replicate the simulation crash at the same density, and why it is happening is a combination of: maybe an over-aggressive timestep (0.001 tau for a dense LJ system with constrained bonds? hmm), and definitely issues minimizing the system with rigid bonds.

During minimization, LAMMPS emulates SHAKE with stiff harmonic bonds – the LJ default is 1e6. This is much larger than the bond strength (100) you were applying. This is not perfect – for it to work, the average bond force (at the applied bond strength) must be much larger than the average total of other forces, and in such a dense system, even one overlap will exert enough force to distort these bonds.

With some quick testing I can see that even at a bond strength of 1e4, the minimization leaves some bonds very far from their constrained length, as seen in this log portion (I’ve attached the input file below):

WARNING: Using fix shake with minimization.
  Substituting constraints with harmonic restraint forces using kbond=1e+04 (src/RIGID/fix_shake.cpp:386)
Neighbor list info ...
  update: every = 1 steps, delay = 0 steps, check = yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 2.122
  ghost atom cutoff = 4.5
  binsize = 1.061, bins = 142 142 1
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair lj/cut, perpetual
      attributes: half, newton on
      pair build: half/bin/newton
      stencil: half/bin/2d
      bin: standard
Setting up cg style minimization ...
  Unit style    : lj
  Current step  : 0
Per MPI rank memory allocation (min/avg/max) = 6.173 | 6.173 | 6.173 Mbytes
   Step          Temp        c_bondmin      c_bondave      c_bondmax        PotEng         KinEng         TotEng         Press
         0   0.013335263    1              1              1              6.1225404e+16  0.0099956597   6.1225404e+16  3.7616888e+16
        10   0.013335263    0.97012697     1.0001073      1.0430835      1.8963e+09     0.0099956597   1.8963e+09     1.1650882e+09
        20   0.013335263    0.84456084     0.99966527     1.207643       1852.7709      0.0099956597   1852.7809      1143.3355
        30   0.013335263    0.83253455     1.0045644      1.5232982      21.112777      0.0099956597   21.122772      4.2179434
        34   0.013335263    0.93753248     1.0035426      1.3783311      6185.6341      0.0099956597   6185.6441     -318.29048
Loop time of 0.0081456 on 4 procs for 34 steps with 2304 atoms

When the dynamics starts, SHAKE must “recover” those bonds, and in a dense system it will likely create new clashes when doing so.

Note that how big a clash will crash the simulation is related to the timestep used – even from the above minimisation, a timestep of 1e-4 (instead of 1e-3) allows for recovery:

   Step          Temp        c_bondmin      c_bondave      c_bondmax        PotEng         KinEng         TotEng         Press
         0   0.013335263    1              1              1              7.1055745      0.0099956597   7.1155702      1.8088773
       500   99.707991      0.999998       1              1.0000018      16.830112      74.737718      91.56783       4.1339008
      1000   43.060012      0.99999821     0.99999999     1.0000017      7.2184397      32.27632       39.494759      1.4855276
      1500   18.807721      0.99999806     1              1.0000017      4.7923341      14.097628      18.889962      0.76162338
      2000   6.7990798      0.99999838     1              1.000002       5.4331775      5.0963589      10.529536      0.24861886
      2500   2.6820699      0.99999822     0.99999998     1.0000016      5.0294338      2.0103883      7.0398221      0.076629885
      3000   1.2987469      0.99999805     1              1.0000016      4.492729       0.97349645     5.4662254     -0.0063997328
      3500   0.63516738     0.99999809     1              1.0000017      4.2535215      0.47609985     4.7296214     -0.017895961
      4000   0.42629755     0.9999984      1              1.0000016      4.04293        0.31953813     4.3624681     -0.016010793
      4500   0.23383998     0.9999981      0.99999997     1.0000017      3.9925577      0.17527849     4.1678362      0.061431975
      5000   0.13920116     0.99999822     0.99999998     1.0000015      3.9576093      0.10434045     4.0619498      0.07905399
      5500   0.063306893    0.99999819     1              1.0000017      3.9552278      0.047452692    4.0026805     -0.0077798396
      6000   0.043964427    0.99999832     0.99999999     1.0000019      3.9334785      0.032954239    3.9664327      0.0956323
      6500   0.027975573    0.99999826     0.99999998     1.0000017      3.9245569      0.020969538    3.9455264      0.042068665
      7000   0.021080426    0.99999808     1              1.0000017      3.9169237      0.01580117     3.9327249      0.0034243914
      7500   0.016066038    0.9999985      1              1.0000018      3.9133556      0.012042556    3.9253981     -0.077301464
      8000   0.014225456    0.99999824     0.99999999     1.0000018      3.9102627      0.010662918    3.9209256      0.021958655
Loop time of 0.622448 on 4 procs for 8000 steps with 2304 atoms

The integrator and equations of motion correctly transfer the high potential energy of the clash into kinetic energy, and the Langevin thermostat successfully dissipates it within 0.8 tau (compatible with the thermostat’s time constant of 0.1 tau). Also note how the eventual equilibrium pressure is not that high – it’s in isolated initial clashes that you are encountering this issue.

Input file:

# presets
dimension      2
units          lj
boundary      p p p
atom_style    bond

# build box
variable scale equal 3 # okay at 1
region         simbox block 0 $(v_scale*50.0)  0 $(v_scale*50.0)  -0.01 0.01
create_box     2 simbox bond/types 1 extra/bond/per/atom 1
mass         *  1.0
molecule        dumbbell mymol.txt
create_atoms   1 random $(v_scale*v_scale*128) 9362 simbox mol dumbbell 8751
velocity        all create 0.01 235092 mom yes rot yes dist gaussian

# potentials
neighbor        1.0 bin
neigh_modify   every 1 delay 0 check yes
bond_style     harmonic
bond_coeff     1 100.0 1.0
pair_style lj/cut 1.122
pair_modify     shift yes
pair_coeff 1 1 1.0 1.0 1.122
pair_coeff 2 2 1.0 1.0 1.122
pair_coeff 1 2 1.0 1.0 1.122
comm_modify     vel yes cutoff 4.5

# thermo and computes
thermo          10
compute bondlen all bond/local dist
compute bondmin all reduce min c_bondlen inputs local
compute bondave all reduce ave c_bondlen inputs local
compute bondmax all reduce max c_bondlen inputs local
thermo_style   custom  step temp c_bondmin c_bondave c_bondmax pe ke etotal press

# fix shake -- only STARTS to give stable results for kbond >= 1e8 if timestep = 1e-3
fix             rigid all shake 0.0001 20 0 b 1 kbond 100

# minimize and run
minimize        1e-4 1e-6 10000 100000
reset_timestep  0
fix             integrate all nve
fix             noise all langevin 0.01 0.01 0.1 543621
fix             two_dim all enforce2d
thermo          500
timestep        1e-3
# does not crash for kbond 1e4 + timestep 1e-4 -- "enough time" to resolve clashes
run             8000
1 Like
10

Dear @akohlmey and @srtee

Many thanks for your valuable comments. Yes, the problem arises from the overlapping of dumbbells using create atoms. I put those dumbbells in a square lattice for my equilibration run. That solves the overlapping problem. Now the production run with shake (without minimization) and input from the equilibration is running fine and I get the delta value 3.7*e-6. That means shake works well now.

11

As some friendly unsolicited advice – please double check your final trajectory with a visualisation. Wouldn’t want that square lattice ordering to persist into your production trajectory and give you skewed results …

12

Thanks @srtee for the suggestion. I will check the production trajectory.