Re: Dielectric constant of water from MSM and PME simulations

From: Zhe Wu (
Date: Fri Oct 23 2015 - 10:46:41 CDT

Hi Mattia,

Just want to put a different issue in a separate email.

Since the eps_RF = infinite may not hold in a 2D system, one could calculate dielectric constant in a more straight forward way: calculate a PMF, the long-ranged part, of two charged particles in your 2D water system. Please find details for such a calculation in my previous paper, page 10528 (last paragraph). <>

But I do encourage you to take a look at the dipole fluctuation formula for dielectric constant calculation or even derive it. I am not sure whether it is difficult to derive but it is not on the top of my head at the moment. I somewhat recall that Charles Brooks might have derived the formula in his 2D-PME papers.

Please keep us updated. It is a very interesting project and I am very eager to learn more!


Zhe Wu
Klaus Schulten's Group,
Postdoc Fellow in NSF's Center for the Physics of Living Cells,
University of Illinois Urbana-Champaign
Beckman Institute, RM 3125,
405 N. Mathews Ave. Urbana, IL 61820
Office: 217-224-3160

> On Oct 23, 2015, at 10:31 AM, Zhe Wu <> wrote:
> Hi Mattia and Axel,
> Thanks for the following up. Thanks Axel for the help! I went through your slides and I really appreciate the very systematic benchmark. Nice work.
> Let’s first focus on Mattia’s first calculation with fully periodic boundary condition. I think the eps_RF = infinite should hold for MSM with FULLY periodic boundary (3D-MSM). So for Mattia’s first simulation, the 3D-MSM result should agree with the PME one.
> The dielectric constant in our JCTC paper was calculated for TIP3P water in CHARMM. There is no question that the dielectric constant of TIP3P-CHARMM is different from the one of TIP3P-EW. In the original TIP3P-EW paper, the different results were shown.
> As Axel suggests, the different results with Mattia’s TIP3P-EW simulations could be just a convergence issue. The reported value of 89 should be converged with the simulation time reported.
> Regarding to Mattia’s second simulation with 2D periodic boundary, I am not sure whether the assumption of eps_RF = infinite still holds. I think Axel could be right and that could be the reason why the dielectric constant calculated with 2D-MSM is so large. Could Mattia try to perform a new MSM simulation with FULLY periodic boundary condition just like the PME one and see what is going on?
> My comments regarding water structure is related to the calculated 258 dielectric constant. I think with 258 dielectric constant, a water molecule could only sense the waters very nearby and thus the water structure could be perturbed. It is effectively the same as performing a water simulation with parameter of eps = 3 instead of eps = 1. Mattia’s g(r) shows a good water structure, which means the water dielectric constant should not that large. The large value could be a result of the dielectric constant calculation only.
> Best,
> Zhe
> Zhe Wu
> Klaus Schulten's Group,
> Postdoc Fellow in NSF's Center for the Physics of Living Cells,
> University of Illinois Urbana-Champaign
> Beckman Institute, RM 3125,
> 405 N. Mathews Ave. Urbana, IL 61820
> Office: 217-224-3160
> <>
>> On Oct 23, 2015, at 9:47 AM, Axel Kohlmeyer < <>> wrote:
>> 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..
>> axel.
>> > 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