From: Otello Maria Roscioni (otellomaria.roscioni_at_unibo.it)
Date: Mon Oct 26 2015 - 12:29:02 CDT
Dear Zhe, dear Axel,
it is very stimulating to receive such a detailed feedback on our
simulations and to discuss them with you.
First of all, I thank Axel for sharing his extensive benchmark of water
properties. Probably the different dielectric constants of bulk water
are a convergence issue. We will extend the production time by 60 ns for
both PME and MSM simulations, which hopefully should yield consistent
results.
Regarding the simulation of a thin-film of water, we wanted first to
exclude that we were doing no obvious mistakes with the simulation
set-up, since it was the first time we used the MSM method. If this is
not the case, we will proceed as follows.
As Zhu suggested, we will simulate a slab of water with MSM plus 3D
boundary conditions. For the sake of simplicity, we will consider from
now on a free-standing film of TIP3P-EW water with two harmonic
potential walls at around 2 nm from the water-vacuum interfaces.
We will also simulate the same system with the PME and MSM-2D methods.
For what concerns the calculation of the dielectric constant, it could
well be that the classic formula for tin-foil boundary conditions does
not hold in the MSM-2D scenario. In this case, I hope the simulations of
the free-standing films with the PME, MSM-2D and MSM-3D methods to
provide some hints.
Meanwhile, we can implement the calculation of the dielectric constant
using the formula 26.3.3D from Prof. Neumann paper for the case where
eps_RF = eps.
Talk back to you soon.
Best wishes,
Otello
On 23/10/15 17:46, Zhe Wu wrote:
> 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).
> http://pubs.acs.org/doi/abs/10.1021/jp1019763
>
> 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!
>
> 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
> http://www.ks.uiuc.edu/~zhewu/
>
>> On Oct 23, 2015, at 10:31 AM, Zhe Wu <zephyrbless_at_gmail.com
>> <mailto:zephyrbless_at_gmail.com>> 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â~@~Ys first focus on Mattiaâ~@~Ys 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â~@~Ys 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â~@~Ys 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â~@~Ys 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â~@~Ys 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
>> http://www.ks.uiuc.edu/~zhewu/
>>
>>> On Oct 23, 2015, at 9:47 AM, Axel Kohlmeyer <akohlmey_at_gmail.com
>>> <mailto:akohlmey_at_gmail.com>> wrote:
>>>
>>>
>>>
>>> On Fri, Oct 23, 2015 at 9:57 AM, Mattia Felice Palermo
>>> <mattiafelice.palerm2_at_unibo.it
>>> <mailto:mattiafelice.palerm2_at_unibo.it>> 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.
>>>
>>>
>>> â~@~Kdear 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
>>> http://klein-group.icms.temple.edu/akohlmey/files/talk-trieste2004-water.pdf
>>>
>>> â~@~Kthe following points matter for computing the static dielectric
>>> constant from dipolar fluctuations (NOTE: not the average but the
>>> fluctuations) :
>>> - â~@~Kit 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, â~@~Ki am not surprisedâ~@~K. 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â~@~Ys 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 akohlmey_at_gmail.com <mailto:akohlmey_at_gmail.com>
>>> http://goo.gl/1wk0
>>> College of Science & Technology, Temple University, Philadelphia PA, USA
>>> International Centre for Theoretical Physics, Trieste. Italy.
>>
>
-- Otello M. Roscioni, PhD Dipartimento di Chimica Industriale "Toso Montanari" Università di Bologna Viale Risorgimento 4 40136 Bologna, Italy
This archive was generated by hypermail 2.1.6 : Thu Dec 31 2015 - 23:22:11 CST