Dear all,
I am experiencing a weird error with lammps. I am trying to relax a box with some large polymers in water. I have circa 7.3 millions atoms. The error says that the exploding bond is between the atoms 7595334 and 7595335. Usually I write a trajectory and every step to investigate the explosion but in this case I can’t and I am not sure if it is a geometry problem or something else.
Could anyone please help me or give me some directions where to look?
I have tried to test different bin sizes and I tried to freeze some parts of the simulations but I always get the same error message at some point.
I am currently using the latest version of lammps. And this is the latest output I obtained:
restoring bond style harmonic from restart
restoring angle style harmonic from restart
restoring dihedral style harmonic from restart
7352445 atoms
4983285 bonds
2824125 angles
874650 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
6 = max # of 1-3 neighbors
12 = max # of 1-4 neighbors
14 = max # of special neighbors
special bonds CPU = 0.181 seconds
read_restart CPU = 1.963 seconds
245070 atoms in group POLYMER
7107375 atoms in group LIQUIDO
35 atoms in group Gruppo1
35 atoms in group Gruppo2
70000 atoms in group Gruppo3
139930 atoms in group Gruppo4
34965 atoms in group Gruppo5
105 atoms in group Gruppo6
2369125 atoms in group Gruppo7
4738250 atoms in group Gruppo8
Resetting global fix info from restart file:
fix style: nvt, fix ID: VeryLong
All restart file global fix info was re-assigned
Neighbor list info …
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 1000000, page size: 100000000
master list distance cutoff = 14
ghost atom cutoff = 14
binsize = 7, bins = 58 58 55
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut/coul/cut, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
Setting up Verlet run …
Unit style : real
Current step : 66840
Time step : 0.5
Per MPI rank memory allocation (min/avg/max) = 1053.0 | 1054.0 | 1054.0 Mbytes
Step Temp Press Volume c_pe1 c_pe2 c_pe3 c_pe4 c_pe5 c_pe6 c_pe7 c_pe8
66840 300.32869 1090.8761 61781284 -564.75453 -150.84135 108030.62 135963.02 -533621.7 131.61584 -20656526 -2501509.8
66841 300.32217 1088.1831 61781284 -589.48253 -151.84304 107968.24 135761.62 -533958.57 139.26338 -20669125 -2499010.1
66842 300.27675 1085.4915 61781284 -619.42795 -144.36151 107690.49 135654.65 -533801.55 139.17748 -20669358 -2504346.8
66843 301.50588 1132.9017 61781284 -622.8919 -145.86747 107648.81 138965.5 -533353.33 133.47996 -20668344 -2496110.7
66844 482.87271 4271.1484 61781284 -639.90236 -153.43334 107431.33 160935.17 -533508.22 132.55837 -20637814 -2494827.1
ERROR on proc 18: Bond atoms 7595334 7595335 missing on proc 18 at step 66845 (…/ntopo_bond_all.cpp:61)
Abort(1) on node 18 (rank 18 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 1) - process 18
And this is the input:
I would say your main problem is that you have not done sufficient testing of your setup on a smaller system, particular of the force field parameters. It is very hard to narrow down a problem when you have millions of atoms. That said, the output you provide (detailed energy info on every timestep) is helpful. I notice that on the penultimate timestep, the temperature and pressure increased by a large factor, clearly indicating an unstable simulation that is scientifically useless from that point on. The error message following this is not necessarily related to the original causation event. Perhaps most significant is the increase in pe_4 by about 20,000 kcal/mol. The simplest explanation is that one of the bonded interactions is too stiff to be resolved by a 0.5 fs timestep. Perhaps you defined it incorrectly. You could try reducing the timestep, but it would be better to do more testing of the force field for a small system. There is probably something unphysical in your energy model.
Aidan
Also, since you apparently have multiple independent instances of this problem, you could use the information contained in :
" Bond atoms 7595334 7595335"
to see if all instances involve the same two atom types in the bond, in which case that would be the smoking gun that would point to a single bad bond coefficient.
Aidan
please note that a statement like “I am using the latest version of LAMMPS” is meaningless for two reasons.
- there are usually 3 different “latest” versions of LAMMPS: the stable version, the regular (or patch) release, and the current development head on github. and those can be different. there can be up to 9 months difference between those.
- what is latest depends on the day. if somebody finds your message later in the mailing list archive because of search for a solution for a similar problem, the term latest has no meaning, because latest will be different at that time.
not to mention that regularly people say “latest” but it is only the latest version available as a precompiled version on their HPC cluster and sometimes those can be rather old.
Dear Axel,
thanks for the note. It am currently using the Lammps stable release from 29 september 2021.
Dear Thompson,
thanks for the help. I tested the force field only on a smaller system only for a short run. I will do a better investigation.
I have a question, why did you think that the bond could be too stiff? What is exactly the problem when a bond is too stiff?
Regards,
Stefano
thanks for the help. I tested the force field only on a smaller system only for a short run. I will do a better investigation.
I have a question, why did you think that the bond could be too stiff? What is exactly the problem when a bond is too stiff?
when doing an MD simulation you do a numerical integration, i.e. you are replacing differentials with finite differences (d becomes delta).
that kind of approximation has limits of how large the delta may be until the approximation breaks down and the finite differences results diverge significantly from the analytical result. you can easily observe this by making a bunch of short test runs with a sequence of increasingly larger timesteps.
how large the timestep may be until the approximation with finite differences breaks down usually depends on how larger the differences are and that depends on how fast particles are moving. what are the fastest moving particles? atoms with the smallest mass and the largest forces. in molecular systems, those are typically hydrogen atoms and the largest forces are from bonded interactions. this is confirmed by experiments: what kind of vibrations do the highest frequencies in infrared or raman spectra correspond to?
people that have seen and done a lot of simulations on molecular systems know from experience that at timestep of 0.5 fs is very much at the limit of what is acceptable for typical molecular force field with unconstrained hydrogen atoms. even if you have simulations that do not have “lost bond atom” crashes, you will have limited energy conservation compared to a timestep of 0.2-0.25fs, which is a conservative choice.
“stiff” bonds, i.e. bonded interactions with large force constants are a particular problem since with those an error in computing the position will lead to much larger forces. just consider a single harmonic bond and you move your atoms in finite increments determined by a fixed time: the faster the motion, the larger the displacement, but that also means that the force is updated less frequently and thus an atom can move higher up the repulsive branch with a larger time step, which will lead to a larger repulsive force, which leads to a larger velocity which leads to larger displacements which leads to moving even higher up the repulsive branch and so on.
there is a lot of research on figuring out ways to have larger timesteps and thus do simulations more effectively and get better sampling of phase space for the same computational efforts. this include schemes like bond/angle constraints (via SHAKE or SHAKE+RATTLE), multi-timestepping (forces for “stiff” interactions like bonds are evaluated more frequently than those for “softer” interactions like non-bonded pairwise forces), mass repartitioning (i.e. mass is shifted from heavier atoms to lighter atoms), advanced langevin style integrators (which involve a friction term that slows down faster moving atoms more).
For more information on the stability of time integration I suggest you read up about the topic in your favorite text book on MD simulations.
Axel.
Dear Axel,
thanks for the answer. I will give a look into that direction.