# Ewald summation without tinfoil bc, but a finite dielectric constant of the surrounding medium

Dear Stan and Chris,

thanks a lot for the helpful reply. Actually, Chris, the system I want to study has a large polarization and Ewald with tinfoil boundaries tries to suppress that. I had a look at “void Ewald::compute(int eflag, int vflag)” in “ewald.cpp” as I guess this is the only method that should be affected by the modifications (forces + energies + virial). Last year I implemented a MD simulation of SPC/E water with Ewald summation. It can be quite hard to spot bugs, that’s why I am reluctant to “touch” the LAMMPS code. However, as you implemented the functionality already I wanted to ask you first, if you would possibly share it, as I have never modified LAMMPS before. If not, I would appreciate, if you could tell me, which methods require modification.

Peter

You shouldn’t have to touch the Ewald routine per se at all, thankfully! All you’re doing is computing the dipole moment of the cell and adding another term to the energy. If you find the bit in the code that does the 2D slab correction it’s very similar to that, just in 3D.

Are you sure about tinfoil boundaries suppressing polarization? Our experience is quite the opposite… the vacuum boundary term acts to counter polarization since the correction term to the energy is always positive and gets bigger with larger system dipole moment. See eg. PRL 103, 207801.

Good luck,
Chris.

Dear Chris,

actually, I think you are right. I just double checked my notes with your paper and the paper from De Leeuw and Perram (1980), Physica 107A (1981) 179-189 and I got the wrong sign for the correction term. For eps_medium < infinity a larger dipole moment results in a higher potential energy, as you said.

Regarding the correction, I am not sure I understood the modifications you suggested. Correcting the energy alone is not enough is it? I mean, U(eps) = U_ewald + f(eps) M^2. Therefore, F(eps)_mC = -grad_mC (U(eps)) = F_mC_ewald - f(eps) * 2* q_mC * M, where mC stands for site C of molecule m, q_mC for the charge of that site, F_mC_ewald is the force you would get with normal Ewald summation and F_mC is the force with the applied boundary correction. To me that looks as if I would also need to correct the force according to the last term -f(eps)2q_mC*M. Am I wrong in this assumption? It would be really helpful if you could tell me where to make these changes, as I still cannot see how it should work without modifying ewald.cpp

Best,
Peter