I have almost finished writing the code of the customary dpd pair style following Axel’s suggestion. Thanks again for Axel’s valuable guidance. But there seems a problem: the computation of conservation force F(i,j) needs the per atom density of both atom i and atom j, but if atom i is in local subdomain and j is not, then even after the reverse communication, the rho(j) is not the correct value. It seems that the reverse communication only updates the local atoms of the subdomain, while if the atom j is in neighboring subdomain, the reverse communication does not update rho(j). It seems more code should be modified to achieve what I want. May be something related to the pack_reverse_comm() and unpack_reverse_comm() function should be modified? What should I do to solve this problem？ Any suggestion is welcome. Since I have no MPI programming experience apologize if this is a naive question.
I have almost finished writing the code of the customary dpd pair style
following Axel's suggestion. Thanks again for Axel’s valuable guidance. But
there seems a problem: the computation of conservation force *F(i,j)*
needs the per atom density of both atom *i* and atom *j*, but if atom *i*
is in local subdomain and *j* is not, then even after the reverse
communication, the *rho(j)* is not the correct value. It seems that the
reverse communication only updates the local atoms of the subdomain, while
if the atom *j* is in neighboring subdomain, the reverse communication
does not update rho(*j*). It seems more code should be modified to
achieve what I want. May be something related to the pack_reverse_comm()
and unpack_reverse_comm() function should be modified? What should I do to
solve this problem？ Any suggestion is welcome. Since I have no MPI
programming experience apologize if this is a naive question.
no MPI programming needed. just look more closely at the EAM code. it
does two communication steps.
the reverse communication that you have already implemented collects data
*from* ghost atoms and sums them to local atoms.
the following forward communication sends data that is stored on local
atoms *to* ghost atoms.
in EAM first the density is computed. when LAMMPS uses newton's 3rd law for
pairs of atoms where one atoms is a ghost and one a local atom, you need
the reverse communication to send that data to the corresponding processor
that owns the ghost atom. this is what happens to "rho". EAM then computes
the "fp" array and that is needed on the ghost atoms as well as on the
local atoms. this is done with comm->forward_comm_pair(this) and the
corresponding pack and unpack routines. if you just need to send the
completed "rho" back to the ghost atoms, you just need to adjust the
pack/unpack for forward communication accordingly.
thus if you follow the EAM example, you should be able to get what you want
without having to write a single MPI command. 
axel.
Thanks, Axel, I did as what you said, and finally get what I want. In retrospect, it is not as difficult as what I thought to make some useful customization to LAMMPS code. And just as what I have heard before, LAMMPS is highly modular and easy to be modified. It is very helpful to me, I hereby express my gratitude to Axel and all developers of LAMMPS.