Equilibrating with Morse Bonds

Cy and Axel,

I have a different theory as to why systems with bond_style morse are
difficult to equilibrate. I've been implementing a 3D version of the
model in [M. L. Falk, Phys. Rev. B 60, 7062–7070 (1999)], where a
number of atomic species interacting by different Lennard-Jones
potentials make up what is similar to a metallic glass. Because I
want to study fracture (and prevent the possibility of cracks
"healing"), I wanted to use a bonding type interaction. The Morse
bond style can be made to look like a L-J potential, and is the only
bond style in LAMMPS where I can create many slightly different L-J
potentials at the same time. Because the systems I'm studying are
amorphous solid, it's not obvious what their density should be. So in
trying to equiilibrate systems with different numbers of atoms per
volume, I ran into the same problems that you appear to be having.

What I realized, though, was that at very low densities, the total
bond energy of the system would approach zero. But the form of the
Morse potential suggests that as (r) goes to infinity, the energy
should converge to (D). If you try the simplest case of two atoms in
a large box separated by various (r), you'll find that the energy at
small (r) decreases to zero at (r0), increases following the Morse
potential above (r0), then abruptly drops back to zero at the
interaction cutoff. So I think it's entirely a non-physical issue
(for lack of a better word), and I'm fairly sure that this leads to
the problems you're seeing.

The easy fix that I came up with was to modify the two lines of code
in bond_morse.cpp that are used to calculate the bond energy. If you
translate the energy by (-D), then the potential naturally converges
to zero at large (r), and this obviously doesn't affect the real form
of the potential in any physically meaningful way. I don't have the
code here to say for sure, but I recall adding "-d0[type]" to each of
these expressions in bond_morse.cpp. I've found that the systems seem
to have much better behavior after that. A nice side effect is that
any bonding will then always lead to negative values of the bond
energy, which is useful in determining what's going on in your system.

I'd actually advocate to the LAMMPS developers to please make this
change in the source code. I could be wrong, but I can't see any
benefit to keeping it like it is.


Thanks for the careful explanation. I think what
you are saying is that the energy function for the
Morse bond has a discontinuity at the pair cutoff
going from D (the formula) to 0.
What I don't understand is
how you can run a system with bonds longer
than the pair cutoff. If the bond atom is not
in the neighbor and ghost list, then LAMMPS won't
find it to compute the bond and should throw
an error. If it does find it (e.g. you run on a single
processor), then it will calculate a huge distance
for the bond, but that should give a energy of
(nearly) D and be fine, since the energy is continuous.

I agree that it is bad for eng minimization
to have a discontinuous function, but I don't
see how that is happening.