Hello,
I am trying to run a simulation of SPC/E water molecules two dimensionally confined between two external planar potentials (Lennard-Jones 9-3 potentials). I would like to add an external electric field to the system along the confining direction so I can calculate the induced polarization os the system, and thereby, the dielectric response.
However, I have the feeling that I am missing something basic that make the calculation wrong. I was reading about the use of the “fix efield” option together with the “kspace_modify slab” option.
I read in another discussion that this is not possible to do with the TIP4P model because fix efield does not properly add the force to this particular model due to the trick to infer the position of the M particle. That is why I am using the SPC/E problem.
However, I also read the following in the LAMMPS webpage: “…If you do it via the fix efield command, it will not give the correct dielectric constant due to the Yeh/Berkowitz (Yeh) correction not being compatible with how works.” If the interaction between the replicas of the system is turned off along the confining direction, why could not we calculate properly the dielectric constant?
My idea is to use the following commands in the input file:
What happens when you do a Kspace calculation with slab correction is the following.
You first compute the regular periodic coulomb interaction in all directions, but then you determine the dipole of the slab system, use Poisson solver to compute the interaction between periodic dipoles in z-direction, and finally subtract that from the regular periodic calculation. That is contrary to your intention to induce a dipole by applying an external electric field. Also, it only cancels dipole-dipole interactions, higher order multipoles still interact across periodic images in z-direction; it is just assumed that those contributions (like dipole-quadrupole or quadrupole-quadrupole) are negligible. But dielectric properties can be “funny” and even the tiniest contributions can make a difference (or not).
I am not an expert in these calculations (@srtee is and may want to chime in), but you may get lucky using MSM instead of PPPM or Ewald. For that you can do non-periodic long-range coulomb without having to apply the correction. Otherwise there may be something useful for you in the ELECTRODE package, but you need to study the docs, examples, associated publications carefully.
Your biggest issue will be in simulating confined water: planar confined water is ordered (up to 9 angstrom from the interface!) and exhibits extremely different dielectric properties from bulk water (see first reference above). Make sure that you are analyzing local properties, and not just the bulk of the system.
Intuitively, I think @akohlmey 's explanation is right. To elaborate on what he has said: the LAMMPS manual recommends slab 3.0 because a certain minimum vacuum distance is needed between images for the dipole correction to work. Presumably, if a larger dipole is induced, this distance will no longer be enough – either to cancel the dipole effect itself, or to cancel the correspondingly larger multipole interactions. Then your simulation cell will exhibit a larger dipole than it would have if it were truly isolated, since it will be polarized by the di/multipoles of neighboring images.
It’s more reliable to use explicit charges, using methods from either of the references in the first paragraph. You could certainly use ELECTRODE to set up those charges (full disclosure: I’m one of the authors and I quite like getting cited), but I don’t think it’s necessary unless you need super-precise measurements (at which point sampling uncertainty and force-field inaccuracy would dominate anyway) or you are interested in fairly non-planar interfaces.
Thank you very much to both of you, the answers were very clear and I really enjoyed reading the suggested papers.
I would like to ask you about another issue which is directly connected with this problem. We had made the same calculations employing the TIP4P/2005 water model. We used the EW3D summations method with fix Efield, substracting by hand the berkowitz correction. The results are surprisingly good if they are compared with other results we have based on DFT. As far as I understood (this was commented in another issue), when fix efield is employed with the TIP4P model, the forces are not correctly applied, instead of applying the force to the massless particle, it is done into the oxygen atom (which has no charge in this particular model), is this correct?
But if this is the case, it is easy to calculate that the torque that would feel each water molecule under the external electric field is 1,35 times higher than the one that would feel the correct TIP4P/2005 dipole with the same external electric field. Thereby, it would be as applying a 1,35 times higher external electrical field (torque = dipole_moment x External field) into the correct system. My big concern with these estimations is that the dynamics of our calculations are not correct, as we are applying the external force into an atom which should not feel it (although it is very close to the massles particle, which is the one that should!). In conclusion, I think that the calculations are not correct, but the results agree surprisingly good, which makes me feel hope :). I would really appreciate to know your opinion about this issue.
Sorry for bothering again, and thank you very much!
To simulate a uniform electric field acting on the virtual charges of TIP4P molecules, you can manually do “field-shifting” by applying one fix efield on the TIP4P oxygens and another on the TIP4P hydrogens, of different field strength, so that the resultant molecular forces and torques are preserved. See pppm_tip4p.cpp:242-252 for some idea of the exact expressions. A third fix efield would have to be applied to non-TIP4P charges.
Legal Disclaimer
This post is intended for educational and informative purposes only and is not intended to serve as mathematical or physical advice. You should consult your supervisor or other local computational chemistry professor before starting this or any other simulation protocol to determine if it is suitable for your unique needs. This is especially true if you or your group have a history of results anxiety or energy non-conservation, or if you have recently experienced Out-of-range atoms - cannot compute PPPM or SHAKE determinant < 0.0 errors when simulating, or have a diminishing quota or grumpy sysadmin which could be made worse by a change in queue submission activity. Do not start this protocol if your advisor or funding administration body advises against it.
It is a bad idea to post a followup question that has little to do with the original question. Please keep in mind that a forum like this also serves as an archive of answers to questions, so that people can find answers that have already been given.
However, to best achieve that, new questions should be posted with a new and descriptive subject line. It is not exactly obvious to expect answers to a question about the TIP4P/2005 water model under a subject that mentions SPC/E.
So please in the future, create a new post with a new subject if the question is not directly related to the discussion topic.
