Simulation Stuck at "WARNING: Shake determinant < 0.0"

LAMMPS LAMMPS Beginners
1

Hello everyone,

I am a pharmacist, studying my Master’s Degree in Medical Biotechnology. I am focused on bioinformatics.

Currently i am working on a project regarding the peptides used in cosmetic preparations and their properties. For this example, i was studying a peptide “Matrixyl” solution and it’s breaking point in thermal stress, and for preparing the simulation i have used the peptide example.

I am sharing with you my data.matrixyl (because of the length kindly access to from my git repo) and the in.peptide file*. I know they are a mass, i just want it to be compiled and get some results (couple of frames for a gif would be perfect).

Also simulation stuck like this:

LAMMPS (3 Mar 2020)
Reading data file ...
  orthogonal box = (16.8402 21.0137 9.7681) to (64.2116 68.3851 57.1395)
  1 by 1 by 1 MPI processor grid
  reading atoms ...
  2004 atoms
  reading velocities ...
  2004 velocities
  scanning bonds ...
  4 = max bonds/atom
  scanning angles ...
  6 = max angles/atom
  scanning dihedrals ...
  15 = max dihedrals/atom
  reading bonds ...
  1365 bonds
  reading angles ...
  786 angles
  reading dihedrals ...
  207 dihedrals
Finding 1-2 1-3 1-4 neighbors ...
  special bond factors lj:   0          0          0         
  special bond factors coul: 0          0          0         
  4 = max # of 1-2 neighbors
  7 = max # of 1-3 neighbors
  15 = max # of 1-4 neighbors
  19 = max # of special neighbors
  special bonds CPU = 0.000424346 secs
  read_data CPU = 0.00683676 secs
System init for delete_atoms ...
PPPM initialization ...
WARNING: System is not charge neutral, net charge = 22.967 (src/kspace.cpp:313)
  using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
  G vector (1/distance) = 0.246118
  grid = 20 20 20
  stencil order = 5
  estimated absolute RMS force accuracy = 0.0319305
  estimated relative force accuracy = 9.61578e-05
  using double precision FFTW3
  3d grid and FFT values/proc = 15625 8000
Neighbor list info ...
  update every 1 steps, delay 10 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 12
  ghost atom cutoff = 12
  binsize = 6, bins = 8 8 8
  2 neighbor lists, perpetual/occasional/extra = 1 1 0
  (1) command delete_atoms, occasional
      attributes: full, newton on
      pair build: full/bin
      stencil: full/bin/3d
      bin: standard
  (2) pair lj/charmm/coul/long, perpetual
      attributes: half, newton on
      pair build: half/bin/newton
      stencil: half/bin/3d/newton
      bin: standard
WARNING: Ignoring 'compress yes' for molecular system (src/delete_atoms.cpp:125)
Deleted 0 atoms, new total = 2004
Finding SHAKE clusters ...
  10 = # of size 2 clusters
  74 = # of size 3 clusters
  0 = # of size 4 clusters
  545 = # of frozen angles
  find clusters CPU = 0.000286066 secs
131 atoms in group peptide
PPPM initialization ...
  using 12-bit tables for long-range coulomb (src/kspace.cpp:332)
  G vector (1/distance) = 0.246118
  grid = 20 20 20
  stencil order = 5
  estimated absolute RMS force accuracy = 0.0319305
  estimated relative force accuracy = 9.61578e-05
  using double precision FFTW3
  3d grid and FFT values/proc = 15625 8000
Neighbor list info ...
  update every 1 steps, delay 5 steps, check yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 12
  ghost atom cutoff = 12
  binsize = 6, bins = 8 8 8
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair lj/charmm/coul/long, perpetual
      attributes: half, newton on
      pair build: half/bin/newton
      stencil: half/bin/3d/newton
      bin: standard
Setting up Verlet run ...
  Unit style    : real
  Current step  : 0
  Time step     : 0.5
WARNING: Bond/angle/dihedral extent > half of periodic box length (src/domain.cpp:906)
SHAKE stats (type/ave/delta) on step 0
  4 7.37924 4.91615 22
  6 7.4382 3.25066 8
  18 1.94236 28.2466 3704
  31 104.015 158.971
WARNING: Shake determinant < 0.0 (src/RIGID/fix_shake.cpp:1742)
WARNING: Shake determinant < 0.0 (src/RIGID/fix_shake.cpp:1742)
WARNING: Shake determinant < 0.0 (src/RIGID/fix_shake.cpp:1742)
WARNING: Shake determinant < 0.0 (src/RIGID/fix_shake.cpp:1742)

Thank you.

2

Apparently not. There is no link here. Also the quoted in.peptide file is not enclosed in triple backquotes (```) and thus it is unusable since the forum software is trying to typeset it.

Please have a closer look at this warning and the SHAKE statistics. Somehow your data file is bogus, since it contains very long bonds. This is confirmed by the fix shake statistics.
For bond type 4 you have an average length of ~7.4 \AA which is nearly 5 \AA longer than the target, for bond type 6 you have ~7.4 \AA with a delta of 3.25 \AA, for bond type 18 a ~30 \AA delta. The warnings about < 0 shake determinant is an indication of degenerate angles which make it impossible to satisfy the constraints.

Please also note this warning. This is very bad. That means either your charge assignments are wrong or that you are missing counter ions to neutralize the system. The net charge should be very close to 0.0.

TL;DR the bond lengths and angles you have defined in your data file are at least in part bogus, plust there is a problem with your charges, which puts the procedure of how you created the data file in question.

3

Thank you for your prompt answer. I have re-edited my post and added my github repo link. Best Regards.

4

Yes, that confirms that the data file is complete garbage. You don’t have proper atom typing, you have bogus charge assignments, the tool creating this file has apparently no knowledge of force fields or the original file that was read didn’t follow necessary requirements.

It cannot even be read with a more recent LAMMPS version due to the inconsistent format.
In the Atoms section the format suddenly changes:

125     1    1    0.00000   17.72349   -1.95412   10.60653 #   C
126     1    2    0.41000   17.73048   -0.98905   10.10059 #   H
127     1    2    0.41000   18.72022   -2.13655   11.01852 #   H
128     1    3    0.00000   17.46071   -2.98310    9.60863 #   N
129     1    2    0.41000   17.51256   -2.56173    8.67911 #   H
130     1    2    0.41000   16.53130   -3.38402    9.69490 #   H
131     1    4   -0.82000   17.27199    5.24865   10.15465 #   O
132     17  14   0.417  57.35121  57.63116  50.79076   0  -1  -1
133     18  13  -0.834  58.09837  59.68005  36.16995  -1   0   0
134     18  14   0.417  58.25901  58.76822  36.41283  -1   0   0
135     18  14   0.417  58.56239  60.19049  36.83355  -1   0   0
136     19  13  -0.834  52.29019  60.51169  50.55611   0  -2   1
137     19  14   0.417  52.61972  60.01708  51.30645   0  -2   1

And when reading this file into VMD with TopoTools (which is more forgiving in terms of image flags) you get something that looks like this:

snapshot
snapshot1045×848 82.4 KB

I suggest you first teach yourself how to do proper force field assignment with either Amber, or CHARMM or OPLS/AA based on available external(!) tools and then convert the resulting topology to LAMMPS. If you don’t have a license for any of those, you can try NAMD and psfgent and the corresponding tutorials, which are available at no cost.

Or - even better - you find yourself a collaborator that has experience in setting up such simulations and and either does it for you or teaches you how to do it. From the files you are providing it is very obvious that you are lacking the knowledge and training to do this on your own.

5

Thank you again, so it is not something to be easily tinkered with by me or someone else. I should start over :slight_smile:
For the tools which you have mentioned, can you suggest any specific courses that i can access?
Best Regards.

6

What tool is best depends on what force field is best and that depends on the details of your system and what you want to learn from your simulations. There is no simple “do this, not that” kind of answer there. Typically this requires some search of the literature and/or consulting with an expert that has worked with similar systems. After settling on a force field, a proper web search should reveal suitable tutorial material. Most of the widely used classical force field packages and their corresponding force fields have been around for a long time (like 20+ years) so there is quite a lot of material around but also a lot of experience. I cannot give a more specific recommendation with the available information.