In my group, we are interested in long range Coulomb forces and a modified Ewald summation technique, so I am studying the long range Coulomb forces. I have seen this post for outputting the long range Coulomb forces, and I have written my own code to do it. I also have my own algorithm that should compute the same forces, but I am getting different results when I look at the autocorrelation functions (same order of magnitude, but different decay behavior). It is certainly possible that my algorithm is faulty, but I want to ask to make sure I am doing things correctly on the Lammps side.
In order to output the long range force, I modified the fieldforce_ik function in the pppm.cpp file, and outputted the qfactor*ekx, y and z on a tagged oxygen particle, which as far as I know should be the long range force on the particle. I also used the gewald command to set the gewald parameter to a value that would be appropriate to compare to my results. I have attached the autocorrelation functions that I computed in this way, which you will see are very different (Ewald split are the forces outputted from LAMMPS). This is a simulation of bulk water in a box of sidelength 3 angstroms, and the tagged particle is an oxygen. Any help in understanding this issue would be greatly appreciated. Thank you!