Compute electric potential when using pppm style


I am trying to compute the long range contribution to the electric potential at various points in a system. To do this I used the compute code from a fellow lammps user (thanks Daniel) modified with corrections for slab geometry (since my system has slab geometry with long range electrostatics (pppm))

Here is how the code does it :

  1. Get the potential values as calculated in the poisson_peratom() function i.e. the u_brick[][][] matrix. ( u_brick*q is the energy per atom )

  2. Extrapolate them to the atom coordinates.

  3. Make the necessary corrections to the term (as done in the poisson_peratom() function skipping the factor q[i] at every step) and that gives the potential at the given point.

Using this compute, I tried running a simple test case, where I place two equal and opposite charges separated by a distance D in the z direction. The system is periodic in the x and y direction.This system mimics a parallel plate capacitor. I use
pair_style coul/long
kspace pppm
kspace_modify slab 3.0
(For now I have kept the cut off small enough such that I do not have any short range contribution to the potential)

Now I place charge neutral atoms at different points inside the box along the z direction. I expect the gradient of the computed potential to be the electric field inside the capacitor, (which would be 4pisigma*C where C is the appropriate factor for the system of units used, and sigma= surface charge density)

As a cross check I also compute the potential energy per atom.
Now for the charged atoms, I get the right value for the calculated potential, and the
pe/atom = q*V
but at the points where I place charge neutral particles the potential is 1/2 the expected value.

(i.e. the field is 2pisigma*C)

I just want to make sure if I am doing something that is conceptually incorrect or is there a problem with the way I set up the system. (using neutral atoms to calculate the potential at various points)

Any help is greatly appreciated.