angle with breakable bonds

Hi LAMMPS users:

I am modifying LAMMPS codes to allow angle potentials to be used with breakable quartic bonds. One thing I need to do is to subtract the pair interaction between two ends of an angle, similar to what is done in the quartic bond potential. However, there are two tricky scenarios while I am coding:

1), an angle shares its two ending atoms with a bond
2), several angles shares the same two ending atoms.

They are possible if one allows one bead to form more than two bonds with other beads. For case 1, subtraction shouldn’t be done in the angle source file. For case 2, subtraction can only be done once. Can someone gives me suggestions on how to make sure subtraction is done correctly without comprising efficiency, because I’ve tried some ways but they really slowed the computation. Or should I even define angles for these two cases?

Best
Shaorui

I think the fundamental problem is to figure
out all the possible angles that the bond is in,
so they call all be "deleted" when the bond breaks.
If you limit yourself to newton on for the bonds,
then the angles are only stored by the central
atom in the angle. So that will have to be one
of the 2 atoms in the bond. The 2 atoms in
the bond could be on different processors, so
you have to make sure they do the same thing,

"Deleting" the angle is somewhat tricky, since
it is stored in permanent and temporary (between
neighbor builds) storage. But the angle potential
presumably knows the 2 end atoms for the
bond that is deleted, so it could compute the
turned-on pairwise interaction.

Hope that helps,
Steve

Hi Steve:

Thanks for the kind help. I’m sorry that I didn’t make my problem quite clear. Here I will explain it with detail.

1). Assume a triangle a-b-c with three bonds a-b, b-c and c-a. In the bond routine the a-c pair interaction is subtracted. In the angle routine I previously wrote, which is incorrect, the a-c pair interaction is subtracted again. Thus a-c pair interaction is subtracted twice, and this should be corrected.

2). Assume a quadrilateral a-b-c-d with four bonds a-b, b-c, c-d, d-a. The a-c pair interaction is also subtracted twice (or even more times if more other angles share a and c as their 1st and 3rd atoms) for angle a-b-c and a-d-c.

To solve problem 1, I tried scanning the bond list, and for problem 2 I tried checking other angles in the anglelist. But this slows the computation a lot (especially for problem 2 since I had to add an inner loop). So is there any smarter way to fix these two problems?

P.S. the variable newton, newton_pair, newton_bond, I know they are Newton’s 3rd law settings and are related to multi-processor computation but still not quite sure about their meanings. Can you please explain it with a little more detail? Thanks.

Shaorui

Hi Steve:

Thanks for the kind help. I'm sorry that I didn't make my problem quite
clear. Here I will explain it with detail.

1). Assume a triangle a-b-c with three bonds a-b, b-c and c-a. In the bond
routine the a-c pair interaction is subtracted. In the angle routine I
previously wrote, which is incorrect, the a-c pair interaction is subtracted
again. Thus a-c pair interaction is subtracted twice, and this should be
corrected.

2). Assume a quadrilateral a-b-c-d with four bonds a-b, b-c, c-d, d-a. The
a-c pair interaction is also subtracted twice (or even more times if more
other angles share a and c as their 1st and 3rd atoms) for angle a-b-c and
a-d-c.

To solve problem 1, I tried scanning the bond list, and for problem 2 I
tried checking other angles in the anglelist. But this slows the computation
a lot (especially for problem 2 since I had to add an inner loop). So is
there any smarter way to fix these two problems?

there may be a way, by doing your own bookkeeping.

you could first create fast lookup table from the bond list
e.g. using the "map". container from stl, or a hash table
implementation (there is one in fix_imd.cpp, that you
could re-use, if you want).

  int b1 = bondlist[i][0] & ~0x0000FFFF;
  int b2 = bondlist[i][1] & ~0x0000FFFF;
  if (b2 < b1) {
   int tmp = b1;
   b1 = b2 ;
   b2 = tmp;
  }
  int key = b1 | b2<<16

  => bondmap[key] = bondlist[i][2];

this way, the smaller atom index number,
will always be in b1 and the key value will
identify this pair of atoms uniquely.
this could handle up to 32k atoms per processor,
which might be sufficient for testing. for production,
you probably want to increase to 64-bit, which
requires a bit more care in programming.

while looping over the angle list, you find out
whether a that is involved in an angle is bond
is broken (type set to 0) by computing keys
for pairs of atoms like above and checking
the bond map/hash for the bond type. that
should be much faster than looping over the
bond list.

at the same time, you compile a second map/hash
for all 1-3 pairs that still have a valid angle interaction
and *not* a bond type >0 or no bond between them
and thus need the LJ term to be subtracted. the unique
keying will avoid double counting.

at the end of looping over all angles, you process
this "active 1-3 angle exclusion" map/hash and
subtract out non-bonded interactions as needed.

the atom index ran be retrieved by reversing the
bitmask/bitshift operations.

P.S. the variable newton, newton_pair, newton_bond, I know they are Newton's
3rd law settings and are related to multi-processor computation but still
not quite sure about their meanings. Can you please explain it with a little
more detail? Thanks.

when doing domain decomposition (and that applies
to a serial calculation in lammps as well), you differentiate
between "local" and "ghost" atoms. "local" atoms are
(more-or-less) inside the domain and "owned" by the
MPI task; ghost atoms are "owned" by the neighboring
domain, or a periodic replica.

the "newton settting" now applies to all atoms where
one is a local and the other is a ghost atom. if newton
is set to "on", then this pair is listed only on one domain
(neighbor and bond/angle/... lists) so forces, energy, and
virial contributions for the local *and* the ghost atom
will have to be stored. with newton set to "off" only the
contribution for the local atom needs to be stored.

is this detailed enough?

axel.

Hi, guys. How is this coming along? Sharoui, have you gotten it working within the current LAMMPS? Are there plans to implement it in the stock version?

Thanks,
Robert S. Hoy
Associate Research Scientist
Departments of Mechanical Engineering & Materials Science, and Physics
Yale University
Phone: (203)-848-7001
Webpage:http://hoymac.eng.yale.edu/~robhoy/