I think a little background will help. Strictly speaking, LAMMPS uses a three-dimensional Ewald sum (or PPPM), which only works for purely periodic boundary conditions in all three directions. There is a 2D Ewald/PPPM method, but it is much more expensive than the 3D version and it is not implemented in LAMMPS. So what is typically done is to use the 3D version of Ewald/PPPM with an additional slab of vacuum, and then correct for the spurious dipole interactions that arise from this approximation (see Daan Frenkel and B. Smit. Understanding molecular simulation from algorithms to applications. Hardcover, November 2002.). So this correction term is the extra term in the energy that you are asking about that is different from the purely periodic case. Thus the “kspace_modify slab 3.0” term you are using is an approximation (you are including only the lowest order correction term and neglecting higher order terms), and can easily be abused if you are not careful. The description in the code is:
Slab-geometry correction term to dampen inter-slab interactions between
periodically repeating slabs. Yields good approximation to 2-D Ewald if
adequate empty space is left between repeating slabs (J. Chem. Phys.
111, 3155). Slabs defined here to be parallel to the xy plane.
Given that warning, the per-atom correction term has to do with the dipole in the z direction, or (qi*zi). So it will be different for the two planes because their position in z is different. In the code, it is:
// per-atom energy
if (eflag_atom) {
double efact = 2.0MY_PIdipole_all/volume;
for (int i = 0; i < nlocal; i++) eatom[i] += qscale * q[i]*x[i][2]*efact;
}
where the term x[i][2] is what makes it different for the two slabs because their position in z is different.
Let me know if you have any other questions.
Regards,
Stan
You can also try kspace_style msm to check your results. It works for slab systems without the use of any correction terms. You should get a similar result as using pppm with kspace_modify slab 3.0.
Stan
Thank you Stan for the detailed description about the slab correction.
However, the results using msm are different from those using pppm in two respects:
-
the values for energy per atom are same for both the z planes in case of msm.
-
the values of energy are different from those reported for pppm (this I can understand could be be because of the accuracy values I choose for pppm and msm as has been the case earlier. )
I can share my scripts if required.
I am trying to include a compute that calculates the electrostatic potential at every atom position in the simulation and was thinking of using pe/atom values as a cross check.
i.e. (short range potential + long range potential) * q = potential energy of the atom (given the only interaction between the atoms is purely electrostatic)
That is why I need an accurate value for the potential energy per atom.
Thanks
Regards
sweta
When you say “different”, can you quantify this, like less than 1%, greater than 30%, etc.? That would be helpful to know if there is a bug in the code.
Thanks,
Stan
So here are the numbers :
(accuracy 1e-4 , rcutoff= 30A, potential energy in kcal/mol )
z plane Using PPPM Using MSM
-50 2146 22.48
150 30.35 22.48
So for z=150, the error is about 35% (assuming the values from msm are correct)
I would suggest shifting your origin so that the two slabs are an equal distance away from the origin, i.e at 100 and -100 in the z direction and try using pppm with the slab correction again, and then compare to msm again. Just a guess, but this may be a requirement for the slab correction that your system violates. Paul can probably comment further. You could also try changing the kspace_modify slab 3.0 to something other than 3.0.
Stan
Yes, that does make the two values equal for the case of PPPM. (since the |x[i][2]| values are the same for the slab correction term for charges in both the planes)
Comparing this to MSM brings the error down to 4. So probably symmetry is a requirement for the way the correction has been implemented or x[i][2] be replaced by x[i][2]-((zhi+zlo)/2)
As for the volfactor value (in the kspace_modify slab), I already checked that one. When the system is not symmetric about the z=0 plane, the difference in pe/atom for the two planes decreases (rather slowly) with an increase in the volfactor. It takes very large values of the volfactor (80-100) for the energies to converge.
Anyway, for now the symmetric system works!
Thanks.
Looks like this is a bug. The force and total energy are translationally invariant, but the per-atom energy is not (my fault–sorry about that). I’ll fix the code and send it to Steve.
Thanks,
Stan