[lammps-users] "Bond atoms ... missing ... " again

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

I am using the 21may2008 version of lammps.
I have some questions about the stability of a MD simulation.
The problem I met is the my MD simulation will stop with error
like "bond atoms missing" sometimes. After checking the mailinglist,
I realize this is a frequently asked "problem". According to the
answers people given, this is not a real problem or bug of the
program itself, but strongly related to the model and parameters
used in the simulation, such as the cutoff, the potential and its
gradient and timestep. I think I have understood these factors.

What I don't understand is how a simulation can become unstable
after a lot of number of MD steps. One of my simulation used CVFF
forcefield and real units. The cutoff is 15 to 20 angstrom, skin size is
2.0 angstrom. The model was converted from Materials Studio. I used
a timestep 0.5 or 0.8 fs, which is normally thought small enough.
According to my experiences, this model should have no problem
to run, and it indeed ran well in Materials Studio and some self-developed
MD proggy. Unluckily, it stopped running in lammps after 300000
MD timesteps, i.e. 150 ps or more. This is unusual to me, because
the properties such as potential energy during the earlier 300000 steps
looks having converged quite well. How can a simulation suddenly go
bad after such long time?

It's true the parameter settings affect the behavior. By decreasing the
timestep or add a 'drag' to the thermostat, the error can be postponed but
never eliminated in my case.

It will highly appreciated if someone can explain this to me.

lin

2

This happens in LAMMPS b/c an atom in a bond cannot
be found by a processor when neighbor lists are built (which
includes lists of bond neighbors). In your case, this
means that the lost atom moved more than 15+2 Angstroms
since the last reneighboring. This could either be b/c
the physics is bad (it should not move that far) or b/c you
are reneighboring far too infrequently.

If the problem is reproducible, I suggest you zoom in on the
timestep that fails (possibly by writing a restart file just before,
then restarting) and doing a lot of output. E.g. thermo output
every step, put a print statement in when reneighboring occurs,
dump the position of the lost atom and force on it every step, etc.
So that you can figure out where that atom has moved to and why.

Steve

3

Some progress of my investigation:

1. the parallel or serial version of lammps produce the same result.
2. The crash invloves an atom which has Lennard-Jones parameters sigma=0.0 and eps=0.0
     This is a 'hn' type hydrogen atom in cvff, which stands for an hydrogen atom bonded to an Nitrogen.
     Automatically assigned by 'msi2lmp'
3. An eye-candy of the atomic details of the crash: http://www.youtube.com/watch?v=d1iWG_4Rp68
    For clarity, many irrelevant molecules are clipped out. Distances between the 'peacebreaker' atoms are
    shown. Two green (a darker and a lighter one) atoms collide then seperate rapidly, but there should have
    no LJ interaction since the darker green atom is a 'hn' with zero sigma and eps.
4. Still need to further investigate the energy calculation and code execution by debugging into the source.
    Still don't understand what is the possible reason and why it happened after a long simulation time.

I know it's hard to get help from others for this specific case. But if someone can give me some hint,
it will be very nice!.

lin

4

If that atom is moving far away from the N it is bonded to, then it
is being acted on by a large force. If its epsilon with everyone
else is 0.0 this shouldn't happen since the only forces
on it will be due to the N bond.

The only exception I can think of is if it overlaps another atom
exactly, and the distance between it and another atom is 0.0.
In that case the pair_lj_class2 code could calculate force = 0.0/0.0
which could be bad. This seems very unlikely since the denominator
would have to exactly 0.0, but you could prevent it from happening
if you set the pair cutoff of that atom (type) interacting with every other
type to 0.0. You could do this with pair_coeff commands or
with neigh_modify exclude commands

I would try to verify that the atom type really has no pairwise interactions
and try to figure out how it is receiving a large force.

Steve

5

Yeah, that's a possible reason. But I have to point out that there is also Coulombic interaction
between the (darker green) 'hn' and the other green(lighter) atom (an O in a c=o group) besides the
'zero' LJ interaction. The hn atom has opposite charge to the Oxygen, meaning they have attractive
interaction. That's probably the reason why they approach each other. All parameters are taken from
the cvff force field. (the hn type there really has an almost zero energy parameter. In the nonbond section
in cvff.frc you will find this line '1.0 1 hn 0.00000001 0.00000'. That's why msi2lmp converts
them into zero sigma and eps in a standard LJ12-6 function form.) This is understandable, because the
hn atom may not need to have LJ interaction with other atoms in order to maintain the correct conformation
and other neighbouring interactions will take care of it. I have no idea if this zero LJ parameters will cause
a problem in lammps.

Another point worthy of an investigation is the Coulombic function, in which an almost zero separation (r)
can be produced when two atoms goes too close to each other. Suppose in the formula q1q2/r, q1 and q2
have opposite sign, and r is extremely small. This will least to huge (overwhelming) attractive force to pull
two atoms going together. However, as said by Steve, two atoms can rarely overlapped exactly, the atoms
will pass by each other, but due to the huge acceleration they may lead to other ill structure such as a over
lengthening bond. Finally, the disaster happened. (Sorry, this is just my imagination, you could take it as
an adventrue movie if it's nonsense to you )

However, I do not this this should happen, because this overlapping problem should be effectively prevented
by the collaborative behaviour of different kinds of interactions (bond, angle etc) if we beleive the force field
is properly translated in a program (including the parameters and function forms).

According to msi2lmp, cvff is defined as Class I ff not Class II. Therefore, the generated data file contains no
cross correlation terms. However, in the original cvff.frc, the parameters fo the cross-terms do exist, such as
bond-bond, bond-angle... msi2lmp doesn't read these. According to the Materials Studio help file, a cvff
has these terms, see http://img13.imageshack.us/my.php?image=ss20090211185941zf3.png .
I rememeber I tried to recover all terms i shown in the picture by manually editing the input data
file and use class2 styles implemented in lammps, but I failed. The reason is the function forms
in lammps cannot meet all terms shown in the picture.In my problematic simulation, I used

bond_style harmonic
improper_style cvff
pair_style lj/cut/coul/cut 15.0
#including 1-4 as nonbonded interaction, same as in MS
special_bonds 0.0 0.0 1.0
#including crossterms, this makes better agreement with Material Studio
angle_style harmonic
dihedral_style harmonic

But to use class2 functions and the cross-term data in cvff.frc, I tried

#with these styles, good energy agreement is found between lammps and MS
#although the angle-angle cross-term is still missing in lammps.
#but the contribution from this term should be small.
bond_style morse
improper_style none
pair_style lj/cut/coul/cut 14.50
#including 1-4 as nonbonded interaction, same as in MS
special_bonds 0.0 0.0 1.0
#including crossterms, this makes better agreement with Material Studio
angle_style class2
dihedral_style class2

With the later settings, for cyclo-Hexane, the energy from lammps is almost the same as calculated in MS using cvff
even I didn't find a function for the angle-angle cross-term (Notice: not the angle-angle-torsion term, which
exists in the current lammps implementation) with a function form E = K * (Theta - Theta0) * (Theta' - Theta0') .
Base on this, I guess my problem may be related to the improper or incomplete use of the cvff force field data.
I breifly check the cvff.frc file, there is actually a non-zero energy parameter under the angle-angle crossterm
'1.0 1 c' n c hn -7.5000'
,in which the hn atom and n and c' are invloved. These are just the atoms involved in my crashed simulation,
although I'm not 100 percent sure it is the reason.

In summary
1. I have to first say I'm not using the 'complete' cvff. This is for sure problematic.
2. If a no-core or very soft-core(softer than the r^-1) could be used with existend of the Coulombic interaction?
Normally the repulsive LJ function will prevent the too-close or overlapping problem ,thus the r in the Coulombic
term will never be too small to lead to a crash. However, as least as seen in CVFF, some atoms can have no
core.
3. How to use cvff completely? Are all the function forms available in lammps? Is there somebody using the complete
cvff and providing some hints?
4. The crash I encountered is rare case, even I used incomplete ff. Just occasionally, the Hydrogen atoms in
a NH-CO group goes close enough to the O atom in the same group, then under the help of attractive Coulombic
interaction they collide and produce huge force due to the small r, although the other angle or bond interaction
prevent them to do so. I have another simulation using the same model and FF, it can run even longer than this crashed one
without a problem.
5. No real debug yet. To know whether the Coulombic interaction causing the problem is easy(almost obvious since
we can see two atoms are almost overlapped in the movie), however, to know what lead them go closer is difficult
because this is a result of many collabrative interactions and also dependent on the local environment.
6. Considering the model can run under other software. Summary 1 and 3 are most probably the reason and the solution.

In this story I must have made mistakes or misunderstandings. Please point them out.Thank you!

lin

Steve Plimpton wrote:

6

I think you need to figure out what 2 atoms overlapped and why
that caused a big force. I don't think it's a LAMMPS problem,
I think it's a problem with your model.

Steve

7

Yes, I know the potentials I used are not fully consistent to what are defined in cvff.
This is the top suspect cause for the crash. I didn't mean lammps is wrong as an
MD integrator, but the function forms implemented for cvff are not complete.
I cannot find a way to include the angle-angle cross-interaction. (if there is, please
point it out for me. )

Without a way to use all the potentials defined in a FF, the simulation can be untable.
One cannot or should not use a part of the force field by simply ignoring some potentials.
I guess everybody here knows this.

I'm sorry if I didn't make things clear. But please give me an answer, could I use all the
potential function forms defined in CVFF?

lin

Steve Plimpton wrote:

8

There is an angle-angle term in the class2 dihedral. There may
be some other term missing somewhere (someone emailed me about this
a while back), but I think this is unlikely to be the source of
your problem. These high order cross terms tend to be very
small for most molecules and should only alter the dynamics slightly.

What is much more likely (in my view) is that the pairwise
interactions are either not set up correctly (by you) to mimic class2
or that LAMMPS is having a problem when 2 atoms overlap strongly,
b/c it doesn't expect the separation distance would ever be 0.0 or
tiny.

However, I don't think your model or the class 2 ff should allow 2 atoms to
overlap, even if LJ epsilon = 0.0 or they have unlike charges, etc. B/c
this is unphysical.

Again, if you can do some debugging (zoom in on the bad timeframe
and dump out more data) to find what 2 atoms overlap and how
the force becomes huge, then I think we could solve this.

Steve

9

Thank you, Steve! I'll debug the problem as you suggested when I have more time.

I think that's me who emailed you about this problem some time ago. But at that time,
I just played a little bit with CVFF in lammps. In the class2 dihedral style, there is a
angle-angle-torsion term, but not the angle-angle term. The angle-angle term is
defined by CVFF (see cvff.frc shipped with lammps tools you will find it.) The difference
between angle-angle and angle-angle-torsion terms is the former one isn't dependent
to the dihedral angle defined by four consecutive atoms. I cannot find a way to use
this function form in lammps without writing a new style. (I hope I'm wrong. Please
tell me how if any one knows. ) Maybe I'm too serious with CVFF. CVFF, who cares?
However, to be honest, there is really no easy way to work with many atomistic
force fields available on the 'market' in lammps, except for the msi2lmp tool. To edit the potential used
for a complex atomistic model manually is painful. Other simulation package like gromacs
has some tools to automatically generate ff parameters for models, which is good. Most of time I work
with coarse-graining models, therefore, anyway, I have to define the inter-bead potentials by
hands, maybe like most of the lammps users doing. This endows large flexibility to lammps,
but in the mean time it also makes lammps not very user-friendly.

I don't know if someone else likes to use a tool like msi2lmp to generate the data input file
as well as the potential parameters. I found this save a lot of work and probably that's why it's
an 'official' lammps tool. However, according to msi2lmp, CVFF is a class I ff, cff91 is class II.
er? Should I use class II styles or other styles in lammps then? In order to get consistent results
with those I got in some other 'prestigious' simulation software, I compared the energy calculation
using CVFF. I found, with carefully selected function forms and some other delicate settings such
as if to include 1-4 interaction, I can exactly reproduce the value in lammps.
(the test is very simple, find another software which can use cvff, and use msi2lmp to generate
lammps input with cvff, then calculate the energy of the model. One should see no difference since
you used the same FF. Is that always true? You have to try it. ) My test result is, for some models, if you are
lucky you could find all the function forms in lammps, you can get the same energy. But most probably
you cannot. The reason is for some models the cross-terms are necessary to get consistent results.
As far as I remember, even you can use some class II styles for cvff you have to manually modify the potential
parameters generated by msi2lmp. (I didn't try cff91. It seems better supported by lammps.)

Sorry for having said some irrelevant words above. I still think lammps is the best.
I'll try to provide some test cases when I'm free.
I first need to debug myself. That's the most difficult mission to accomplish for a human being. :slight_smile:
lin

Steve Plimpton wrote:

10

I don't know a lot about class 2, but I think the LAMMPS implementation
is missing at most one cross term. If you zero that one out, then you
should be able to get the energy of a LAMMPS system and an MSI system
to agree. I know people have done that test.

Ditto for other force fields. I agree the conversion into LAMMPS format
can be tedious or painful. But the tools provided should work where they're
applicable. For example, the ch2lmp tool takes CHARMM output and converts
it into LAMMPS input and people do the energy test you describe with good
results. The msi2lmp tool may be somewhat out-of-date - I don't know many
people who use it, and I don't maintain it, since I didn't write it.

Steve

11

Thanks, Steve! I have to verify the force field before I can use them anyway.
I can also implement the missing function form if necessary. So I will give some feedback once I have the results.

lin

Steve Plimpton wrote: