Question on fix shake internals


I’m currently trying to borrow an idea from fix_shake for my own fix development. Basically, I need to do something like what the angle version of shake does, and get data about a small group of atoms which are connected together, do some calculations with that information, and then update the data based on the result of these calculations.

What fix_shake appears to do in the angle case is to use atom->map() on a set of known tags for the atoms of interest, and then uses the resulting ids from that function to find distances and angles on all processors (since there isn’t any selection criterion that I can see). Then, once the calculation is done, it updates the correct forces only on the proc where the ID falls within its nlocal set (the if(i0<nlocal) clause at the end, for example).

My problem is that when I try to do something similar, I get segfaults, which I suspect is the result of trying to use invalid IDs from atom->map() on one of the processors. In fact, it seems like what fix_shake does should segfault as written, since it appears to use the returned IDs on all processors, even though based on the force update criteria later, it looks like those IDs should be out of the array range for processors that don’t actually own the atoms.

I was hoping someone could explain what’s going on here; either what I’m apparently misunderstanding about the angle function in fix_shake, and/or how exactly atom->map() works and how to correctly use and interpret its return values in a case like this.

Thomas Allen

All atom->map() is convert a global atom ID into a local
index on that proc. It can return -1 if the proc doesn’t

know that atom (which would be an error in fix shake).

The local index can be for an owned atom (< nlocal) or

ghost. The latter is OK for fix shake to update, e.g modify
a force on, since those forces on ghost atoms can be
communicated back to the owning atoms (later in the