What is the difference between the per-atom array in atom class and the fix->array_atom?

Dear Axel,

I am coding a custom fix to compute a custom per-atom array in every N timesteps. I checked fix_store_forces.cpp and fix_property_atom.cpp, and it seems that these two fixes are using different mechanisms in LAMMPS:

In fix_store_forces.cpp, the per_atom_flag is set to 1, size_peratom_cols is set to 3, and per_atom_freq is set to 1. The above three flags are inherited from the class fix and the usage is self-explanatory. Then the array_atom pointer (also inherited from the class fix) is assigned by each local atom. I think I can understand this and this is exactly how I did in my own fix, and it worked well.

However, in fix_property_atom.cpp, the non-standard per-atom array is created through atom->add_custom() and then atom->add_callback() is called. I checked the code for atom->add_callback() but it is beyond my knowledge.

Here are my questions:

  1. In the first method, I would like to ensure the ghost atom information is correctly communicated. What is the minimal set of virtual function in the class fix that I need to override? Currently I override copy_arrays(), pack_exchange(), unpack_exchange().

  2. What is the different between fix::pack_exchange() and fix::pack_comm()? I checked the documentation and the code comment but the explanation is somehow vague to me.

  3. In the second method, what is the usage by calling atom->add_callback()? I tried adding this code in my custom fix, and something strange happened: the same code can run on MacOS, MPICH-4.2, gcc-11, 10 cpus, but having segmentation fault error on Linux, same setting, 80 cpus. I tried to use gdb to debug but it only tracks down to srealloc(), not very useful.

  4. What is the advantage of fix->array_atom over atom->add_custom()? I believe the developer designed this for a reason. Just curious.

With all respect, your reply is greatly appreciated.

The main purpose of fix property/atom is to implement per-atom properties that can be references with either i_name or d_name and so on in a number of places. This is a very special process with specific requirements. It should not be used for regular per-atom data implemented and managed by a general fix.

You send per-atom data from local atoms to ghost atoms, you need to do a so-called “forward communication”. This is implemented by overriding the functions “pack_forward_comm()” and “unpack_forward_comm()” and then calling “comm->forward(this)” when you want to send the data.
As an example you can look at the fix_group.cpp source code, which needs to update the atom->mask information on ghost atoms after a dynamic group update.

There is no Fix::pack_comm() only Fix::pack_forward_comm() and Fix::pack_reverse_comm().

Fix::pack_exchange() is for an “exchange communication” which migrates atoms and their associated data between processors. This is usually done before rebuilding the neighbor lists and communicating ghost atoms.

The Fix::pack_forward/reverse_comm() are used for forward communication (send per-atom data from local atoms to ghost atoms) and reverse communication (sum data added to ghost atoms to their corresponding local atom. Needed when newton_pair or newton_bond are enabled).

This function needs to be called (typically with the Atom::GROW argument) so the fix is registered with the Atom class to allocate and initialize additional data in case new atoms are created with the system. This is rarely needed for fixes that compute per-atom data. They can usually grow their per-atom data as needed when the fix is called. The callback is only needed if the per-atom data needs to be always present. Most fixes, especially those that update their data every step, call “memory->grow” on their own.

Don’t use it! It is not meant for your use case.

See my previous explanation. They serve different purposes. Fix property/atom can “augment” an atom style to provide a per-atom property not normally included with this atom style in use or custom data that can be used in many ways. It is a method to avoid having to program a new atom style.

Thank you for the very rapid response! The explanations are mostly clear to me.

One following question: Now I see pack_exchange() and unpack_exchange() is not for communicating ghost atoms. Could you give me a concrete example when they will be called? I assume one case is, when the per-atom array maintained by my custom fix is set as output by dump custom as f_FixID, the above two function will be called because the MPI rank 0 will receive all data from every other processor then make the output.

Please let me know if my understanding is correct or problematic. Thank you! And have a good night!

Your assumption is wrong. Only “owned”, i.e. local atoms will be communicated. Most fixes that compute properties, only do so for local atoms. Ghost atoms are only needed when the fix data is used internally, e.g. from a pair style. Thus all ghost atom related communication would be unnecessary for them.

Exchange communication is required for fixes that have “persistent” data, i.e. information that may be used at any time and has to be also present when the fix computation is called and thus it has to persist when local atoms migrate between subdomains. If the fix’es per-atom information is always computed from scratch it is not needed.

Border communication is needed to create and update ghost atom data from local data.

Forward communication is needed to update ghost atoms with newly computed local data.

Reverse communication is needed to collect data stored with ghost atoms to their corresponding local atoms.

Thanks for correcting me and this is very helpful.

Exchange communication is required for fixes that have “persistent” data, i.e. information that may be used at any time and has to be also present when the fix computation is called and thus it has to persist when local atoms migrate between subdomains. If the fix’es per-atom information is always computed from scratch it is not needed.

I am still confused with this “persistent” concept. Since the fix per-atom data is managed by the dynamic pointer Fix->array_atom, isn’t the life cycle of this variable span all the timesteps that invoke this fix, and then available at anytime or anywhere by using this pointer?

No. The data is only valid on time steps where the fix was invoked. Otherwise atoms may have migrated and the fix data is would be outdated and the array sizes incorrect. Take fix store/force for example. Its data is “recomputed” in every step, so it does not need to be “persistent” and per-atom data migrated. This is in contrast to fix ave/atom where the data is updated only on every “Nfreq” timesteps, but has to be accessible on all timesteps so further data can be averaged. So the exchange communication support is needed in case atoms migrate between subdomains in between updates of the averages. However, it does not require data on ghost atoms.

The data is only valid on time steps where the fix was invoked.

I checked the code of fix store/force and fix/ave/atom. Indeed, if a fix updates its data in every timestep, comm is not necessary.

I assume the word “valid” means array_atom is in correct size (at least larger than atom->nmax) and in correct order (array_atom[i] corresponds to the local atom i). So if I do not override Fix::pack_exchange() and Fix::unpack_exchange(), and my fix only update the array_atom by Nfreq timestep, the data will probably be wrong during other timesteps. Is this correct?

Based on the pseudocode of LAMMPS timesteps, I guess Fix::pack_exchange() and Fix::unpack_exchange() will be called after the reconstruction of neighbor list, at the comm->exchange() step. Please let me know if this is right or wrong.

Thank you!

Yes. LAMMPS should be able to detect whether the fix is updated on a compatible time step from the value of Fix::peratom_freq. This is, for example, done by dump style “custom”, which will throw an error, if the fix wasn’t recomputed on the time step of the dump.

Wrong. It is called before building the neighbor lists. Here is the relevant block of pseudocode.

nflag = neighbor->decide()
  if nflag:

First it is determined whether a neighbor list needs to be build at all (Neighbor::decide()). If not, only a forward communication is needed to update data on ghost atoms from the local atoms.

Second a bunch of preparation steps are done. Some fix hooks run, ghost atom data deleted, box and sub-domain info updated and the neighbor list binning (sorting of atoms into bins to speed up the build) prepared.

Then, local atoms that have moved outside their subdomains are exchanged with the corresponding neighboring subdomain and the list of ghost atoms recreated.

Only now, the neighbor lists can be rebuilt (Neighbor::build()), followed by some post-neighbor list update hooks for fixes.

Thanks Axel, for the detailed explanations. It helps a lot.

Hi Axel,

Another question: I found that local atom sorting can cause data mismatch in my custom fix per-atom data. That is: I update Fix::array_atom by every Nfreq steps, and output the array by the fix custom command, every Noutput steps. If I leave the local atom sorting on, data mismatch will happen if Nfreq is bigger than Noutput. If I turn the atom sorting off by atom_modify sort 0 10, this mismatch disappears. So I assume this is due to my code does not deal with the atom sorting correctly. How can I make it behave correctly like fix_ave_atom.cpp? I already override the pack_exchange() and unpack_exchange() functions but it seems that they are used for communications between processors, but my problem here comes from local processor internal sorting.

Thank you!

Well, the obvious answer to that is: study what fix ave/atom does.

I checked the fix ave/atom code. Interestingly the only difference is fix_ave_atom.cpp called atom->add_callback(Atom::GROW) in its initiator. I guess this is the way that fix need to use to tell atom class to rearrange fix’s per-atom array.

I originally added atom->add_callback(Atom::GROW) in my code, and then it somehow caused error on HPC, as I mentioned. Now I moved this function call to the very end of the fix initiator and it magically works on both HPC and my laptop. Have no idea for why but the problem seems solved. The data mismatch problem disappeared.