Bond quartic question

Dear List,

I'm somewhat puzzled about the implementation of bond quartic.

The documentation:
> If a bond length ever becomes > Rc, LAMMPS “breaks” the bond, which means two things.
> First, the bond potential is turned off by setting its type to 0, and is no longer computed.
> Second, a pairwise interaction between the two atoms is turned on, since they are no
> longer bonded. LAMMPS does the second task via a computational sleight-of-hand.
> It subtracts the pairwise interaction as part of the bond computation. When the bond breaks,
> the subtraction stops. For this to work, the pairwise interaction must always be computed by
> the /pair_style/ <> command, whether the bond is broken or not. This means that /special_bonds/ <>
> must be set to 1,1,1, as indicated as a restriction below.

The code does exactly this. But for two bonded beads, why is LAMMPS calculating a WCA force contribution in the bond, then subtracting a WCA pair force again via force->pair->single(..) in the bond, and finally adding it again via the pair calculation?

With special bonds 111, it would be much simpler just to define the bond potential WITHOUT a WCA contribution, and just switch this off the bond when it is broken, and hence always calculate the WCA term as part of the pair interaction independently of whether beads are bonded or not.

The only difference, I can see, is that the partitioning of energies between Epair and Ebond would differ when bonds are broken, but this is a philosophical question. But it would make more sense with Ebond =0 in the limit where all bonds are broken.

Not sure I understand. I don’t think there is a requirement that the LJ term in the

bond style formula be the same as the pairiwise LJ defined for the 2 bonded atoms.

The user is free to set the sigma,epsilon for the latter to whatever they wish.

In fact it’s not required that the pair style be LJ, it could be anything. I’m not saying

all options make sense, but they are allowed. Which implies doing the addition/subtraction?

Also, the bookkeeping of epair vs ebond might be imporant for some analyses?

The current bookkeeping is correct I think? The bond formula is tallied with ebond,

the pairwise subtraction is untallied from epair. And ebond -> 0.0 when all bonds

are broken.


Dear Steve & List,

Agreed that the complexity is needed if bonded WCA differs from non-bonded WCA.

It's because I'm interested in breaking cross-links between Kremer-Grest polymers, and using bond/quartic seemed overly complicated and computationally expensive in this case. I'm implemented a simplified version, which assumes special_bonds 111 and bonds are just FENE springs without the WCA contribution.