Slab issue

Hello,
I am doing simulations of SPC/E water between two hydrophobic (LJ 9/3) walls with an electric field (.5 V/A) applied to the water molecules in the z-direction. For the kspace commands, I use the commands

kspace_style pppm 1.0e-4
kspace_modify slab 10.0

and have the correct boundary conditions (fixed in the z direction). I am using an nvt thermostat with 100.0 fs. relaxation time.
I have been calculating the polarization potential for these simulations

Polarization Potential (z) = -1/(epsilon_0) \int_{-L}^{z} \int_{-L}^{z’} dz’ dz’’ \rho^q(z’’),

where the boundary is at -L. According to electrostatics, the slope of the polarization potential should satisfy the equation

\epsilon = e_field/(e_field-slope)

because the water reduces the applied field by a factor of \epsilon. When I calculate the value of \epsilon this way, I am getting a factor of about 100, but it should be around 70 for SPC/E water.

People in my research group have done this calculation using DL POLY and gotten the correct results in the past. I may be implementing this simulation incorrectly, but I have been trying to find bugs in my code for about a month, and am starting to think there might be a subtle issue with the implementation of slab corrected 3d Ewald summation techniques from the kspac_modify slab command.
Has anyone else faced issues such as these. Do you have any ideas as to what might be going on, or are you aware of any known bugs that might be causing this error. I would really appreciate any assistance! Thank you.
-Teddy Baker

teddy,

there are not any known bugs related to your setup. i have not done
the kind of calculation you are doing, so i cannot comment on that
specifically.

what is a bit unusual are the choice of settings. kspace_modify slab
10.0 is very large. the poisson solver applied requires the vacuum to
be at least twice the size of the charge distribution, so using a
value of 3.0 is already ample and allows for some diffusion. since you
have confinement walls, there should be no concern about that,
provided the walls are correctly places. on the other hand, using pppm
1.0e-4 is a rather sloppy convergence parameter. so i wonder, if
choosing it more tightly, has some impact on your calculartion.

axel.

Hi again,
I have done a number of simulations more to test how various parameters affect the results I am getting, and unfortunately the result has remained stubbornly the same. I have changed the accuracy to very tight values (10e-8), have tried using full Ewald summation, have changed the Coulomb cutoff radius, and have seen the same results for the dielectric constant each time.
It seems that I might be doing the analysis incorrectly, but I have double checked my code a number of times, which amounts to computing the charge density through a histogram and integrating it to find the average electric field in the bulk. Any more suggestions for things to try would be greatly appreciated!
-Teddy

Hi again,
Also, I forgot to mention that I also tried changing the vol_factor for the slab correction, which also did not improve the result.
Just to recap the problem, when I calculate the value of the dielectric constant for SPC/E water confined between walls I get a value that is different than that in the published literature by about 30 (100 vs. 70). I am not sure if this is a problem with my analysis, my setup, or the slab geometry Ewald algorithm.
-Teddy

Stan may wish to comment on this thread …

Steve

Hi again,
I have run a number of simulations, and the problem still persists. I wrote my own algorithm for the Yeh and Berkowitz slab correction, and in fact the results look better, although they still give a value for the dielectric constant that is a little bit higher than the correct value. I tried this simulation with much better accuracy (10e-8), using ewald summation, and using pppm ewald. The best results I have gotten so far are with my algorithm and full ewald summation, although the simulation runs very slowly.
I have attached some plots showing the slope of the polarization potential in the bulk, which is the average electric field produced by the water in response to the external field. You can see that there is a marked difference between my algorithm (labeled YB for Yeh Berkowitz) and the Lammps implementation, and the best is my algorithm using full ewald summation (some simulations were have not completed yet).
Also, when I ran the simulation with the Lammps slab algorithm for a smaller box, it ran extremely slowly, I am somewhat puzzled by that.
I am starting to think that there is indeed a bug in the Lammps implementation of the slab geometry.
-Teddy

Slope_Compare.png

Slope_Small.png

Hi Teddy,

What is different between your algorithm for the Yeh and Berkowitz slab correction vs the one in LAMMPS? It is a simple correction, so it should be easy to tell. Your graphs are very noisy–are you sure the differences are statistically significant?

Are your results sensitive to the Lennard-Jones cutoff or system size effects? Are you using similar values for LJ cutoff and system size vs simulations you are comparing to in literature? Have you checked energy conservation/timestep size? Many things could affect your results…

Stan

Hi Stan,
For my algorithm, I added the dipole correction of Yeh and Berkowitz and manually changed the box size to a larger length, while keeping the walls the same (had to disable an error message in fix_wall.cpp to do that). I looked at your code, and it seems you add the same dipole correction, so my guess is there is something different in the way the Ewald sum is implemented.
I’m sure the results are statistically significant, I have ran these simulations many times and get a slope that is noticeably off, yielding a result for the dielectric constant that is consistently wrong (~100 instead of 70). I have tried using a larger LJ and Coulomb cutoff, that did not fix the problem. I am using a timestep of 1 fs., the results in the literature I refer to use the same timestep. I did not check energy conservation.
If there is another thing I haven’t checked I am more than happy to do so, I have been trying to solve this bug for over two months and it is quite an obstacle for my research, so any other suggestions would be greatly appreciated!
-teddy

Try calculating the dielectric constant using fluctuations without a field–I have done it before for a different system in LAMMPS and it isn’t hard. See if you get the expected value. Then turn your field back on and calculate the dielectric constant using fluctuations again. See if it gives you the expected value or if it matches your “wrong” results from your other method. Then you can narrow down the problem. Here is a reference: http://dx.doi.org/10.1016/j.cplett.2011.02.064

Hi again,
I think that that approach might be a little bit round about given that there is already published literature that has done that, so I already know the correct value, which is not the result I’m getting.
-Teddy

But you don’t know why you aren’t getting the correct value. Did you implement your SPC/E water model correctly in LAMMPS–does it give you the correct value for a bulk periodic system without a field, walls, and slab correction? You won’t know for sure until you test it. We did something sort of similar using the slab correction in LAMMPS, and it worked fine, see Figure 13: http://dx.doi.org/10.1103/PhysRevE.91.022309.

Stan

That is a good point, I may be implementing the SPC/E water or the simulation incorrectly. However, it will take me a long time to calculate the dielectric constant from an equilibrium simulation, and if I get the wrong result I will not have pinpointed the problem in my simulation.
Instead, do you have an example code of SPC/E water that you could send me to see if I am doing everything correctly, or alternatively if I attached my simulation code could you check it and tell me if you see anything wrong? That would help me a lot.
-Teddy

Well, there is this: http://lammps.sandia.gov/bench/bench_spce.tar.gz . I didn’t implement the model, so I can’t say if it is correct or not.

However, it will take me a long time to calculate the dielectric constant from an equilibrium simulation, and if I get the wrong result I will not have pinpointed the problem in my simulation.

If you got the wrong results, you would know that it didn’t have to do with the slab correction and was something more fundamental. That seems like pinpointing the problem to me. J

Stan

   That is a good point, I may be implementing the SPC/E water or the
simulation incorrectly. However, it will take me a long time to calculate
the dielectric constant from an equilibrium simulation, and if I get the

i don't think that statement holds. the static DC is not system size dependent.
you can compute it from tiny bulk systems containing as little as 64
water molecules.
please see pages 18 to 21 of the following PDF:

http://klein-group.icms.temple.edu/akohlmey/files/talk-trieste2004-water.pdf

axel.

Hi again,
I looked at the example code you sent me, which I think I actually based my spc/e simulation on (in addition to a different example). Two things I noticed

1.) The charge I use is .4238 for H but there it uses .4236 (I got my value from here http://www1.lsbu.ac.uk/water/water_models.html). I’m not sure if that would make a big difference, but as far as I can tell I am using the accepted value.
2.) I do not append the nx ny nz information at the end of the atom coordinates for the read in file, I’m not sure if that would make a difference.

Other than that, everything else is exactly the same. Do you think item 2.) would make a difference?

I can check the dielectric using the equilibrium code, however if I got the correct result I think that would prove that there is a bug in the Lammps slab implementation as far as I’m concerned.
-Teddy

Oh I think that probably I don’t need the nx ny nz information because I am not using the replicate command, but tell me if that is incorrect.

Like I said, given that my entire file is essentially the same as the one you sent me, I am inclined to think that my simulation is not the issue.

Ok I’m going to compute the dielectric constant using the bulk equilibrium simulation and my code. If the result is correct then I think that is strong evidence that the slab implementation has an issue.

I have verified the ewald slab implementation by comparing to an analytic solution. Just place two oppositely charged particles in a slab system. They form two infinite sheets of charge due to periodic boundary conditions in the x and y directions. If the charges are far apart enough in the z-direction, the force in the z-direction converges to 2piq^2/(Lx*Ly), (see Understanding Molecular Simulation, Frenkel and Smit, p. 318), and LAMMPS reproduces this result.

Stan

Hi again,
I think that if there is a bug it would be a fairly small one, because the calculation of the dielectric constant is a fairly subtle calculation (it depends on the cancellation of two large values to produce a relatively small value). I think most of the results you would get are correct, this is probing the dynamics very sensitively.
I’ll let you know how my equilibrium results turn out, the simulation is running pretty fast so I should know in a few days.
-Teddy

Hi all,
I’ve gotten about 5 ns. of data for the bulk water calculation of the dielectric constant (it took me a few runs to debug my code for calculating the mean square dipole moment). I’m getting a value for epsilon of about 73.4, which I think is within acceptable margin of error (given the slow convergence you showed me in the powerpoint).
This seems to suggest that my raw simulations are not the problem. It could still be that I’m implementing the slab correction or wall incorrectly (I don’t think the electric field is the problem), or that my data analysis is wrong. I have checked my data analysis a number of times though, but of course something could have slipped by.
-Teddy