[lammps-users] Need efficient method for summing densities over all neighbors within force cutoff

I am developing pair_eam_poly.cpp for calculating EAM potentials using a polynomial h(x) as in the PRL paper by Alfredo Caro. To calculate elemental concentrations, I need the individual densities of each element, summed over nearby atoms, for each atom i. I originally tried to do this inside the loop (copied from compute() in pair_eam.cpp) which computes densities at each atom, but several of the densities rho were zero. Evidently the way the Spatial-Decomposition method is implemented does not list neighbors all the way out to the cutoff distance. So I made a separate loop to sum densities at each atom in the simulation, over all neighbors in the simulation. This achieves the desired result, but at the expense of slowing the execution significantly. Is there a better way to do this? I.e. is there a slick way to sum densities of all neighbors j within the cutoff distance of atom i, without exhaustively visiting every atom in the simulation twice?

     -- Al


Have you tried your first approach with the newton
flag set to off? See:


With newton set to "on" (the default), some
interactions are communicated from neighboring
processors rather than explicitly computed. So that is
probably why you were seeing the zero densities in
some cases.

There is probably a more elegant solution to this
problem, but switching newton off may provide a quick


The normal EAM computes the total density at each atom, using
the 6-way stencil communication routines in LAMMPS. This
is done correctly for newton pair on or off. You should be
able to do something similar for an array of quantities (e.g. density
for each atom type).