[gmx-users] viscosity calculation using g_energy (3.3.3)

[JMandumpal]

 Dear Berk and David,

                                          Thank you very much for your appropriate and informative replies.  I tried another method (traverse current method) to calculate the shear viscosity ( a non equilibrium method, which has been described in Berkś paper : Journal of Chemical Physics, 116, page 209 ( Determining the shear viscosity of model liquids from molecular dynamics simulations)),   

                           I used the g_tcaf utility (ie g_tcaf -f traj1.trr -s binary.tpr -oc test.xvg) . As suggested by David, I increased the system size ( from 500 to 2048 TIP4P molecules). I ran in NVT ensemble which allows the pressure to fluctuate.
Apart from that I added following options to my mdp file, where accelaration of 1A/ps² was given to the system.

;NON EQUILIBRIUM STUFF
acc_grps                  = system
accelerate                = 0.1 0.0 0.0
cos_acceleration          = 0.1

----------------------------------------------------------------

Moreover, I saved the trajectory in every 1ps ( so total 500 frames for a 500ps simulation)

then,

I got the following output:
k  1.593  tau  1.000  eta  0.09835 10^-3 kg/(m s)
k  1.593  tau  1.000  eta  0.09835 10^-3 kg/(m s)
k  1.593  tau  1.000  eta  0.09835 10^-3 kg/(m s)
k  2.252  tau  1.000  eta  0.04917 10^-3 kg/(m s)
k  2.252  tau  1.000  eta  0.04917 10^-3 kg/(m s)
k  2.252  tau  1.000  eta  0.04917 10^-3 kg/(m s)
k  2.252  tau  1.000  eta  0.04917 10^-3 kg/(m s)
k  2.252  tau  1.000  eta  0.04917 10^-3 kg/(m s)
k  2.252  tau  1.000  eta  0.04917 10^-3 kg/(m s)
k  2.759  tau  1.000  eta  0.03278 10^-3 kg/(m s)
k  2.759  tau  1.000  eta  0.03278 10^-3 kg/(m s)
k  2.759  tau  1.000  eta  0.03278 10^-3 kg/(m s)
k  2.759  tau  1.000  eta  0.03278 10^-3 kg/(m s)
k  3.185  tau  1.000  eta  0.02459 10^-3 kg/(m s)
k  3.185  tau  1.000  eta  0.02459 10^-3 kg/(m s)
k  3.185  tau  1.000  eta  0.02459 10^-3 kg/(m s)

---------------------------------------------------------------------

Which shows a strong k dependence over the property: shorter k, better the viscosity, as pointed out in the paper. However, the value obtained is around 0.01 times less than the experimental value (1pa-second). Adding to that, the results obtained by this method seems to be very convincing unlike the g_energy that shows a great divergence!!

So the situation is getting better now. Now, I would like to know whether this can be improved if I save the trajectories more frequently ( 500 fs) and run for longer, say 2ns or change value of accelaration . 

Any thoughts ?


regards,
Jes.


On Thu, 12 Mar 2009 Berk Hess wrote :
>
>Hi,
>
>This is a very inefficient method for determining the viscosity.
>Also you need really perfect pressure fluctuations: NVT, shifted potentials,
>probably even double precision.
>There was a mail about this recently.
>There are better methods, have a look at:
>http://dx.doi.org/10.1063/1.1421362
>
>Berk
>
>Date: Thu, 12 Mar 2009 07:39:52 +0000
> From: jesbman at rediffmail.com
>To: gmx-users at gromacs.org
>Subject: Re: Re: [gmx-users] viscosity calculation using g_energy (3.3.3)
>CC:
>
>
>David,
>
>
>
>Thanks for the quick reply.
>
>
>
>Indeed I did  as what you suggested- g_energy -f water.edr -vis test.xvg
>
>
>
>The output file created includes three columns.
>
>
>
>1. time ( ps) 2. shear viscosity (3) I assume it is bulk viscosity.
>
>
>
>It seems, the unit given is cp. ( 1cp= 1* 10¯3 Pascal Second).
>
>
>
>The bulk viscosity of water at 300 K is approximately 0.7 cp. But the value ( Bulk viscosity) I got from the program gives me 100 pa-s, an increase of two order of magnitude. I wonder whether I have done anything wrong while specifying the frequency of saving energy file.
>
>
>
>I have saved the energy file in every 2ps. Isn´t that enough for a simple system like water? OR should I have to save trajectories in every 5fs as suggested by one in a previous post.
>
>
>
>I post the first 20 lines of the output file.
>
>
>
>-------------------------------------------------------------------
>
>
>
># This file was created Thu Mar 12 16:20:09 2009
>
># by the following command:
>
># g_energy -f water.edr -vis test.xvg
>
>#
>
># g_energy is part of G R O M A C S:
>
>#
>
># GROup of MAchos and Cynical Suckers
>
>#
>
>@    title "Bulk Viscosity"
>
>@    xaxis  label "Time (ps)"
>
>@    yaxis  label "\8h\4 (cp)"
>
>@TYPE xy
>
>@ view 0.15, 0.15, 0.75, 0.85
>
>@ legend on
>
>@ legend box on
>
>@ legend loctype view
>
>@ legend 0.78, 0.8
>
>@ legend length 2
>
>@ s0 legend "Shear"
>
>@ s1 legend "Bulk"
>
>    1.99203      9.6633     96.3893
>
>    3.98406     11.1625     98.1365
>
>     5.9761     12.6631      99.838
>
>    7.96813     13.4652     101.366
>
>    9.96016     13.7012     100.249
>
>-------------------------------------------------------------------------
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090313/bde1997e/attachment.html>