Re: Dielectric constant of water from MSM and PME simulations

From: Axel Kohlmeyer (
Date: Fri Oct 23 2015 - 09:47:38 CDT

On Fri, Oct 23, 2015 at 9:57 AM, Mattia Felice Palermo <> wrote:

> Dear Zhe,
> thank you very much for your accurate answer.
> > Could you specify how you calculate the dielectric constant?
> It has been computed with a inhouse program using the classical
> expression for the dielectric constant from the average dipole moment
> <M^2> as in your JCTC paper.
> > The dielectric constant of TIP3P-EW water is 89 (D. J. Price, C. L.
> > Brooks, J. Chem. Phys. 2004, 121, 10096), while you are getting 96.8
> > with PME and 83.3 with MSM. You can find the way we calculated the
> > dielectric constant on page 773 of our JCTC paper.
> We suspect that the higher value found in our simulation with PME might
> be related to the size of the sample (N=1000 in the TIP3P-EW original
> paper against N=11000 in our simulations). However the value is not far
> from the literature value of 89 (TIP3P-EW, model B) and from the
> value of 104 reported in your JCTC paper. We also checked other
> properties e.g. the oxygen-oxygen radial distribution function and
> the average module of the molecular dipole moment, and they are in
> agreement with the expected results for the model.

​dear mattia,

some comments on this subject. it has been a long time, since i worked on
dielectric properties of water as a graduate student and later got a chance
to do some follow-up work in my spare time.
for your reference, please have a look at pages 18-21 of

​the following points matter for computing the static dielectric constant
from dipolar fluctuations (NOTE: not the average but the fluctuations) :
- ​it is a *local* property. system size has very little to no impact. when
you sit down derive the expression of the fluctuation formula, you will
quickly see that the number of molecules will cancel out.
- it is a *slowly* converging property. most published data is simply not
fully converged.
- it is of monumental importance to use the correct fluctuation formula.
most likely you have used the one for ewald summation with conducting
boundary condition. in that case what you are simulating is effectively an
infinitely large sphere assembled from rectangular unit cells embedded in a
conducting dielectric (epsilon = infinity). i seriously doubt that this is
applicable to MSM. it would be more likely, that you need to derive the
fluctuation formula for a situation where epsilon is equal to the computed
epsilon. but don't take my word for it. these issues are extremely subtle
and i had to twist my brain quite hard for months to finally get a grip on
it and come up with consistent and convincing explanations for the
observations from my simulations.

on the notion of g(r) being reproduced. that means *nothing*. you need to
get them right, too. but with the same g(r) you can have very different
total system dipole fluctuations. the conducting boundary conditions of
commonly used ewald summation (and equivalent) actually enhances
fluctuations. if you switch to a vacuum embedding (epsilon = 0), then those
total dipole fluctuations are strongly suppressed, yet you can't tell a
difference from the g(r), as the average structure is practically
unaffected (you are looking at the second moment, after all, and the
average cancels out).

> > It is very surprising to see a dielectric constant of water to be
> > 258. How is the water structure? With such a large dielectric
> > constant, the water structure should be very unphysical.
> The first surprising fact is that the dielectric constant of bulk
> water differs if the MD simulations is carried out with PME and MSM
> methods. The dielectric constant of water with value 258 was obtained
> for a thin film of water supported on silicon dioxide and using the
> MSM method, while when using the PME method we obtain a dielectric constant
> of 91.5. Visual inspection of the samples through VMD did not show any
> unphysical structure of the water. We also computed the oxygen-oxygen
> radial distribution g(r)_O-O for the films, which you can find attached,
> and
> it does not show any significant difference between the two simulations.

based on my observations outlined above, ​i am not surprised​. and more
importantly, you cannot easily transfer the method to other simulation
environments. you will need to derive the proper fluctuation formula for
such a situation as well..


> > Also, could you give us a little bit more details about your VDW
> > cutoff scheme? As reported in TIP3P-EW’s original paper, a model
> > that incorporates a long-range correction for truncated VDW will
> > give a dielectric constant of 76 instead of 89. Are you using the
> > same cutoff scheme for both PME and MSM?
> Yes we used the same cutoff scheme for both simulations. More in detail,
> these are the settings for the MSM simulations:
> cutoff 12.0
> switching on
> switchdist 10.0
> pairlistdist 14.0
> timestep 2.0
> nonbondedFreq 1
> fullElectFrequency 2
> stepspercycle 10
> langevin on
> langevinDamping 1
> langevinTemp 298.15
> LangevinPiston on
> LangevinPistonTarget 1.01325
> LangevinPistonPeriod 250
> LangevinPistonDecay 100
> LangevinPistonTemp 298.15
> cellBasisVector1 63.340 0.0 0.0
> cellBasisVector2 0.0 63.340 0.0
> cellBasisVector3 0.0 0.0 63.340
> rigidBonds all
> MSM yes
> The same setting were used for the PME simulations, except of course for
> the calculation of electrostatics:
> PME yes
> PMEGridSpacing 1.5
> Thank you very much for your help!
> P.s.: I changed the subject of this email to have all future
> contributions sorted in the same thread. Sorry for the confusion!
> --
> Mattia Felice Palermo - Ph.D.
> Università di Bologna
> Dipartimento di Chimica Industriale "Toso Montanari"

Dr. Axel Kohlmeyer
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.

This archive was generated by hypermail 2.1.6 : Thu Dec 31 2015 - 23:22:10 CST