Problems with dihedral_style potential

LAMMPS LAMMPS Mailing List Mirror
1

Dear Community,

I am a Master student from Germany and new to MD simulations with LAMMPS. My current goal is to simulate a coarse-grained polymer system, i.e. a system of M chains each made of N beads using the attached force-field.

Unfortunately, the LAMMPS simulations are always crashing after a few minutes or, sometimes, hours (ERROR on proc 0: Bond atoms 366 367 missing on proc 0 at step 2865170 (../ntopo_bond_all.cpp:63)).

The problem seems to be related to the dihedral_style potential, because when I am not using any (or instead use dihedral_style zero) the simulation finishes normally.

Now, I understand that if the bonds in a dihedral become collinear, the corresponding dihedral angle can not be defined any more, which should cause significant numerical errors in the simulation. Or in my case I checked certain output step by step using:

"compute 1 all property/local aatom1 aatom2 aatom3",
"compute 1 all property/local datom1 datom2 datom3 datom4", and
"compute 1 all property/local patom1 patom2" (with all the appropriate computes and dumps)

and

a) there are no exact 180° or 0° angles (only occasionally 179.99)

b) there are no noticeably strange dihedral values (such as NAN)

c) there is always at least one exceptionally high Potential (Lennard-Jones) Energy value between two beads of different chains (note that I am using up to 3'rd neighbor exclusions) in the last times-step before the crash.

Because of this I would exclude the above mentioned reason for my crashes.

So it seems there is a correlation betweeen the dihedral_style potential and the overlapping of non-bonded beads. But how can there be a causality between these two events and how can I prevent it?

Cheers

Simon

ff_short.txt (655 Bytes)

2

Dear Simon,

The error says that bond number 366 and atom number 367 is missing. It must be probably due to following reasons:

  1. Polymer segment around particular bond must have overlapped due to strong dihedral forces. (Indicates poor choice of dihedral coefficients and cut-offs).
  2. Your simulation box isn’t maintaining its dimensions (due to high volume fluctuations) and since your boundaries are periodic, as your box side reduced to dimension less than your polymer unit’s dimension, overlapping may occur along that particular dimension. However this must not be the case with your system, since you aren’t applying any typical thermodynamic constraints yet, as per your input file.

However, I would strongly recommend not to set dihedral interactions to zero because that will cause significant error in your final desired output, since polymers do undergo significant scissoring, rocking and other modes of vibration. Thus dihedral interactions are significant in majority of non-planer units.
As per me, I would suggest you to optimise your choice of dihedral coefficients and cut-offs, since I feel there must be problem 1 as mentioned above.

Deep

3

The problem seems to be related to the dihedral_style potential, because
when I am not using any (or instead use dihedral_style zero) the
simulation finishes normally.

Actually, I think Simon was already going in the right direction here. See my comments below.

Now, I understand that if the bonds in a dihedral become collinear, the
corresponding dihedral angle can not be defined any more, which should
cause significant numerical errors in the simulation. Or in my case I
checked certain output step by step using:

“compute 1 all property/local aatom1 aatom2 aatom3”,
“compute 1 all property/local datom1 datom2 datom3 datom4”, and
“compute 1 all property/local patom1 patom2” (with all the appropriate
computes and dumps)

and

a) there are no exact 180° or 0° angles (only occasionally 179.99)

You just discovered the source of your error.

It’s not necessary for θ = 180 in order for a numerical explosion to occur. Numerical explosions occur whenever the force acting on a particle is large compared to what the timestep you are using can support. (This in turn is usually chosen to be considerably smaller than the period of oscillation of the bonds or bond-angles you are using.)

The magnitude of the force due to dihedral interactions (which is proportional to ∂φ/∂x) diverges as θ -> 180. (If I’m not mistaken, it diverges proportional to 1/sin(θ), where “φ” is a 4-body dihedral angle, and “θ” is one of the associated 3-body bond angles, and “x” is one of the coordinates of the atoms, the force on which you are calculating.) So if θ = 179.99, (and if I didn’t make a mistake), then the dihedral force would have been increased by a factor of approximately 5729.6 (≈180/(0.01π)), compared to what they would have been if θ=90 degrees.

I suggest using dihedral-potentials which have a strong barrier within 20 or 30 degrees of 180.0. It can be useful to use “dihedral_style table” for this:

https://lammps.sandia.gov/doc/dihedral_table.html

Comment: But if you are getting θ=179.99, then I suspect you are using a bond-angle force with an equilibrium bond angle θ0=180 degrees. If that is the case, then there is no barrier preventing θ->180. In that case, it simply does not make sense to apply a dihedral angle force to this polymer. I suspect you might need to increase the complexity of the representation of the polymer. A simple chain of beads is often not a good model of a polymer which can support twisting. I suggest reading the simulation papers comming from Andrzej Stasiak’s lab (Try reading the appendix from this paper.)

On a related note, I created a dihedral-style which makes it possible to avoid the singularity as θ -> 180. You can read about it here:

https://lammps.sandia.gov/doc/dihedral_spherical.html
I use this dihedral style in several polymer models, including this one (image1, image2):

b) there are no noticeably strange dihedral values (such as NAN)

179.99 is very strange. (see above)
In summary: Unless you are using dihedral_style spherical (see above), it is very inefficient, and probably dangerously unstable to allow θ to ever approach within 20 degrees of 180.

c) there is always at least one exceptionally high Potential
(Lennard-Jones) Energy value between two beads of different chains (note
that I am using up to 3’rd neighbor exclusions) in the last times-step
before the crash.

That would be a result of these atoms moving too close together, likely because of the large force that the dihedral angle exerted on one of them, causing it to move suddenly and overlap with another atom. I suspect when this happens, you probably encountered a θ value near 180 degrees a few timesteps beforehand.

Good luck.

Andrew

4

Hello Andrew and Deep,

This problem also shows up, when I am lowering my timestep to less than a femtosecond… Okay, I understand this. I already read about this dihedral_style, but so far I am not able to reproduce my polymer characteristics with that…

5

Comment: But if you are getting θ=179.99, then I suspect you are using a bond-angle force
with an equilibrium bond angle θ0=180 degrees. If that is the case, then there is no barrier
preventing θ->180. In that case, it simply does not make sense to apply a dihedral angle force
to this polymer. I suspect you might need to increase the complexity of the representation
of the polymer. A simple chain of beads is often not a good model of a polymer which can
support twisting. I suggest reading the simulation papers comming from
Andrzej Stasiak’s lab (Try reading the appendix from this paper.)

Yes, that’s true! Two of my angular potential have an equilibrium at 180. However, the exact
same simulation (with tabulated potentials) was done in GROMACS before.
There, it worked perfectly.

It sounds as if the GROMACS implementation of your model did not have this problem, but the LAMMPS conversion does. In that case, you probably made a mistake when you converted it into LAMMPS format. Mistakes like that are really easy to do.

All I can say is if any of the 4-body dihedral angle (φ) interactions in the model contain 3 atoms which (θ) approach 180°, then you are guaranteed to eventually have a numerical explosion in your simulation, regardless of the timestep. (Incidentally, 1fs is not that small, at least for all-atom simulations.)

(It is possible to have a twistable polymer containing 180° bond-angles. As an example, the coarse-grained twistable polymer model that I mentioned earlier from Andrzej Stasiak’s lab has some 3-body bond-angles constrained to 180°, however those 3 atoms do not appear in any of the dihedral interactions. [2 of the 3 atoms do, but not all 3.] This is difficult to explain without drawing pictures, but I’m too lazy to do that. See the appendix of that paper.)

I had the task to export the Simulation to LAMMPS to check if
we could get a simulation speed up.

Compared to Gromacs? Seems unlikely. (I don’t know. LAMMPS has other nice qualities: flexibility, generalizability, numerical accuracy perhaps.)

First, I copied the tabulated potentials from GROMACS into LAMMPS, but with
LAMMPS the simulations used to crash after some time. Then, I changed from
tabulated potentials to LAMMPS potentials, but that did not help.
Does someone have an idea, why this problem is not showing up with GROMACS?

Perhaps when you converted the model into LAMMPS format, some of your angle or dihedral interactions are connected to the wrong atoms. Again this kind of mistake is easy to do. (Also check the two ends of the polymer to see if anything weird is happening there.)

I suggest building a small version of your system.
Try a polymer containing only 2-4 monomers and containing the minimal number of dihedral interactions (hopefully only 1).
See if the problem persists. If not, then keep adding monomers until the problem returns.

In addition, I would draw a cartoon stick-figure of your short polymer on paper and see where all of the 3-body angle and 4-body dihedral interactions are located. Make sure that none of the 3-body angle interactions with 180 equilibrium values contain 3 atoms which are participating in the same 4-body dihedral angle.

If the GROMACS model (or the script to convert it into LAMMPS) was written by somebody else, then for the sake of your own understanding, perhaps you should try building the system from scratch for LAMMPS instead of relying on somebody else’s work. (Assuming you are not already doing so.) It should not be that hard. The LAMMPS DATA file format is relatively easy to understand. You can write a script to create these kinds of files relatively easily.

Alternatively, you can try using one of the LAMMPS molecule-builders to do this for you. The topotools molecule builder might be useful to you, if you like VMD and don’t mind using TCL. I will also include a shameless plug for my own tool which I wrote (when I was in your position) to build general coarse-grained molecules and polymers in LAMMPS. (Examples: 1, 2, 3.) If your polymer is complex enough, or if you eventually plan to make the system complicated and add other molecules, then in the long run I suspect it might be less cumbersome to use this tool than to write your own script. You can use these instructions (1, 2, 3) to make the polymer wrap to any shape you desire (and then add additional bonds or branches in moltemplate if you need them). You can use moltemplate without sacrificing any understanding about what your system is or how it behaves. You’re still in control. It is not a black box tool. (Neither is topotools.)

(It actually seems that https://www.sciencedirect.com/science/article/pii/S0032386119304707?dgcid=rss_sd_all (Kubo et al., 2019),
had the same problem using LAMMPS)

If so, it is not because of a bug in LAMMPS.

Good luck!

Andrew