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.