When you say poisson_peratom() in your step 3, I am not sure what this means (this function doesn’t have any factor of q[i]).
I would suggest using a particle with charge=1 to get get V as a cross-check from pe/atom = q*V, instead of using neutral particles. BTW for two point charges in a box, the per-atom energy for a single atom is only half the total interaction energy of the 2 interacting particles, which may be where your factor of 1/2 comes from, not sure though.
Stan
Sorry Stan, I should have been more careful.
Its not poisson_peratom(), I include the corrections as calculated in the compute() function in pppm.cpp where corrections to the energy per atom are calculated.
For example this term in the original code :
g_ewaldq[i]q[i]/MY_PIS + MY_PI2*q[i]*qsum /
(g_ewaldg_ewaldvolume);
is modified to :
g_ewaldq[i]/MY_PIS + MY_PI2qsum /
(g_ewaldg_ewaldvolume);
in my compute that calculates the electric potential.
Yes, you are right, the per atom energy is half the total interaction energy and that is what I observe for the point charges at either end of the box. So the potential for all charged atoms in the box is correct and agrees with the pe/atom values.
If I add an additional charge=1, to get the V at any point in this system, , the system would no longer remain charge neutral, that is the reason I add neutral atoms.
Since the results for charged atoms is correct, i was wandering if for neutral atoms, I am over correcting the calculated potential or if there is something I am missing in the calculations.
Here are my results for the simple system that I use:
(Energies in Kcal/mol and potential in Kcal/(mol e) where e represents charge.
Total Energy of the system : 25.03
q x y z Potential pe/atom expected value of potential
-1 0 0 0 -12.5156 12.5156 -12.5156
0 0 0 48 -0.534976 0 -1.0700
0 0 0 49 -0.267487 0 -0.5397
0 0 0 50 1.00075e-15 0 0
0 0 0 51 0.267487 0 0.5397
0 0 0 52 0.534976 0 1.0700
1 0 0 100 12.5156 12.5156 12.5156
Thank you
Regards
sweta
I’d start with a periodic system (instead of slab). Add two point charges fixed in space, say q1 = -0.7 and q2 = -0.3, and then add a third charge of q3 = 1.0 and move it around. Then you can compare your computed potential with the modified pppm code to the per-atom potential for charge 3 (for a system that is charge neutral) and make sure that they agree, or see if you are off by a factor of 2 again. Then try changing the values of the charges, try it for slab, etc.
Stan
I tried the different cases with and without slab. As long as there is some charge on every particle in the system, the potential and the pe/atom values are in agreement, for both the periodic system and slab geometry.
So I am assuming there is no problem with the way I calculate the potential or at least that is what the results suggest.
Now coming back to the parallel plate capacitor example, I tried the following :
I replaced the q1=-1 at (0,0,0) and q2=+1 at (0,0,100) charges with q1=-1, and q2=+0.99 respectively. And now I move around a third charge q3=+0.01 along the z axis (instead of having neutral atoms)
Since the test atom is now charged, I cross checked it with the value of pe/atom and pe/atom=q*V.
However when using the computed potential to calculate the electric field inside the capacitor, the field is still off by a factor of 2. Somehow only one of the charges q1 or q2 is making a contribution to the field at q3. I am not sure what am I missing out!
Results as shown :
3 charges with slab
Total Energy of the system : 24.82749 Kcal/mol
q x y z Potential pe/atom expected potential value
-1 0 0 0 -12.4607 12.4607
0.01 0 0 50 0.069800 0.00069800 0.13368
0.99 0 0 100 12.4911 12.3661
Total Energy of the system : 24.8327 Kcal/mol
q x y z Potential pe/atom
-1 0 0 0 -12.4624 12.4624
0.01 0 0 51 0.33545 0.0033545 0.66576
0.99 0 0 100 12.492 12.367
Regards
Sweta
How are you calculating the potential from the parallel plate capacitor? Shouldn’t it just be a linear function of distance and voltage from the two plates? For example, if I take your original system of only 2 point charges, you said that you expected the potential to be -0.5397 for a z position of 49:
Total Energy of the system : 25.03
q x y z Potential pe/atom expected value of potential
-1 0 0 0 -12.5156 12.5156 -12.5156
0 0 0 48 -0.534976 0 -1.0700
0 0 0 49 -0.267487 0 -0.5397
0 0 0 50 1.00075e-15 0 0
0 0 0 51 0.267487 0 0.5397
0 0 0 52 0.534976 0 1.0700
1 0 0 100 12.5156 12.5156 12.5156
But if I use a linear combination of -12.5156 + (49/100)(12.5156 - -12.5156) = -0.25031, which is close to what the code gave (not off by a factor of 2). Is that not right?
Stan
Yes you are right, if I use the results from the simulation as the basis that makes sense,
but I use the electric field to calculate the potential at a point.
So for the case of a parallel plate capacitor, the field due to plate 1 and plate 2 is :
E1=2pisigma1 / epsilonr
E2=2pisigma2 / epsilonr
where sigma = surface charge density, epsilonr is the dielectric constant
So the potential at a z plane,
V = (E1* d1 + E2* d2)
In the example, sigma1 = -0.01 e/A^2, sigma2 =0.01 e/A^2, epsilonr=78, d1 = 49-0 =49, d2=49-100=-51
V= ((2pi0.01)/78)*(49-51) e/A
= -0.0016 e/A
in ( kcal/mol e ) units this would be -0.53497.
And that is the value I would have expected the potential to be. Did I miss something here? I must have, because if the value of the potential at either end is correct (which it is ), then the slope of the potential has to be 0.26… and not 0.53…but I can’t see where the difference is coming from.
Thanks
Regards
sweta
Take a look at http://physics.stackexchange.com/questions/65191/what-is-the-electric-field-in-a-parallel-plate-capacitor and http://faculty.wwu.edu/vawter/PhysicsNet/Topics/Capacitors/ParallCap.html. I think you may be double counting the field.
Stan
No, I think the field calculation is correct. According to the second link : E (in the center) = 4pisigma =(0.534 kcal/mol e)
With reference to this example
3 charges with slab
Total Energy of the system : 24.82749 Kcal/mol
q x y z Potential pe/atom expected potential value
-1 0 0 0 -12.4607 12.4607
0.01 0 0 50 0.069800 0.00069800 0.13368
0.99 0 0 100 12.4911 12.3661
Total Energy of the system : 24.8327 Kcal/mol
q x y z Potential pe/atom
-1 0 0 0 -12.4624 12.4624
0.01 0 0 51 0.33545 0.0033545 0.66576
0.99 0 0 100 12.492 12.367
Essentially what I did was moved a charge of 0.01 by 1 Angstrom.
So using the change in energy of the two systems (dEnergy) , the value of charge (q3) and the change in z (dZ) the electric field can be calculated as :
dV = dEnergy / q3 = (24.8327 - 24.82749) / 0.01 = 0.5209
Electric field = - dV / dZ = -0.5209
Another method that gives the expected electric field value.
Not sure about the field, but I don’t think there is a bug in the code. For example, if you sum up the per-atom kspace energy, it equals the total kspace energy. How to interpret the results is another matter. Your question is more a physics question than a LAMMPS question… 
Stan
BTW, it looks like you are off by a factor of 2 everywhere, including at the edges of the box, using your previous example. Say you look at at z=100:
V = (E1* d1 + E2* d2)
In the example, sigma1 = -0.01 e/A^2, sigma2 =0.01 e/A^2, epsilonr=78, d1 = 100 =49, d2=0
V= ((2pi0.01)/78)*(100) e/A
= -0.0016 e/A
in ( kcal/mol e ) units this would be -26.748, not -12.5156.
So you probably need to multiply all your potential values by a factor of 2. This probably has to do with the way that the per-atom energy is defined: as half the interaction energy between a pair of interacting particles.
Stan