[gmx-users] Problem with g_dipoles and SPCE water

[Tom Kirchner]

On 04/10/2012 12:00 PM, gmx-users-request at gromacs.org wrote:
> On 2012-04-10 11:16, Tom Kirchner wrote:
>> >  On 04/10/2012 11:10 AM,gmx-users-request at gromacs.org  wrote:
>>> >>  On 2012-04-10 09:57, David van der Spoel wrote:
>>>>> >>>  >    On 2012-04-10 09:43, Tom Kirchner wrote:
>>>>>>> >>>>  >>    Hi all,
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    I want to test the g_dipoles tool and made a 20 ns simulation with SPC/E
>>>>>>> >>>>  >>    water. Sadly the result for the dielectric constant (e_sim=109.713)
>>>>>>> >>>>  >>    looks very different from experimental (e_exp=78) and reported values
>>>>>>> >>>>  >>    for SPCE (e_paper=70). I put the command line and the output below.
>>>>> >>>  >
>>>>> >>>  >    Did you run the simulation with PME or at the very least RF?
>>>>> >>>  >
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    Gromacs 4.5.3, SPC/E water, 20ns, NPT, 300K, 1bar
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    g_dipoles -corr total -c dipoles -f md_20ns.xtc -s md_20ns.tpr -b 0 -e
>>>>>>> >>>>  >>    20000
>>>>> >>>  >    You don't need the corr flag for computing the dielectric constant.
>>>>> >>>  >
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    #==========================================
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    Using 5 as mu_max and -1 as the dipole moment.
>>>>>>> >>>>  >>    WARNING: EpsilonRF = 0.0, this really means EpsilonRF = infinity
>>>>>>> >>>>  >>    Selected 0: 'System'
>>>>>>> >>>>  >>    There are 1024 molecules in the selection
>>>>>>> >>>>  >>    Using Volume from topology: 30.2049 nm3
>>>>>>> >>>>  >>    Average volume over run is 30.2103
>>>>>>> >>>>  >>    t0 0, t 20000, teller 10001
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    Dipole moment (Debye)
>>>>>>> >>>>  >>    ---------------------
>>>>>>> >>>>  >>    Average = 2.4883 Std. Dev. = 0.1342 Error = 0.0000
>>> >>  Another point: the literature value of the dielectric constant for SPC/E
>>> >>  is for the rigid model. Looks like your using a flexible model here.
>>> >>
>>> >>
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    The following averages for the complete trajectory have been calculated:
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    Total<   M_x>   = 1.45986 Debye
>>>>>>> >>>>  >>    Total<   M_y>   = 5.71202 Debye
>>>>>>> >>>>  >>    Total<   M_z>   = 3.39196 Debye
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    Total<   M_x2>   = 10964.8 Debye2
>>>>>>> >>>>  >>    Total<   M_y2>   = 10794.4 Debye2
>>>>>>> >>>>  >>    Total<   M_z2>   = 10762.4 Debye2
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    Total<   |M|^2>   = 32521.6 Debye2
>>>>>>> >>>>  >>    Total |<   M>|^2 = 46.2638 Debye2
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    <   |M|^2>   - |<   M>|^2 = 32475.4 Debye2
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    Finite system Kirkwood g factor G_k = 5.12208
>>>>>>> >>>>  >>    Infinite system Kirkwood g factor g_k = 3.43028
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    Epsilon = 109.713
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    #==========================================
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    I would be really grateful if anybody could point me to my error. I
>>>>>>> >>>>  >>    attached my .mdp file in case the error lies there.
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>    Best regards
>>>>>>> >>>>  >>    Tom
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>
>>>>>>> >>>>  >>
>>>>> >>>  >
>>>>> >>>  >
>>> >>  -- David van der Spoel, Ph.D., Professor of Biology
>> >  Thank you for your reply,
>> >
>> >  I used the rigid model. I also used PME, I think this is the error, for
>> >  g_dipoles seems to be constructed for reaction field. Is this correct?
>> >  I'll redo the simulation with reaction field now.
>> >
>> >  Best
>> >  Tom Kirchner
>> >
>> >
> No! g_dipoles does not care. You may have a remnant command in your mdp
> file (-DFLEXIBLE), please remove it. For  rigid model there should not
> be any fluctations in the dipole.
>
Many Many thanks for your help, there was an error in the topology file.

I am wondering, why the flexibility of SPC/E water has such a big effect 
on the dielectric constant. Also, which formula or method does g_dipoles 
use to calculate the dielectric constant. In the mailing list you refer 
to a JCP paper of 1998 by you, but there only a formula for reaction 
field is given.

Best
Tom
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20120412/58f870c4/attachment.html>