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.