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