<P>
Dear Berk and David,<BR>
<BR>
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&#347; paper : Journal of Chemical Physics, 116, page 209 ( Determining the shear viscosity of model liquids from molecular dynamics simulations)), <BR>
<BR>
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.<BR>
Apart from that I added following options to my mdp file, where accelaration of 1A/ps² was given to the system.<BR>
<BR>
;NON EQUILIBRIUM STUFF<BR>
acc_grps = system<BR>
accelerate = 0.1 0.0 0.0<BR>
cos_acceleration = 0.1<BR>
<BR>
----------------------------------------------------------------<BR>
<BR>
Moreover, I saved the trajectory in every 1ps ( so total 500 frames for a 500ps simulation)<BR>
<BR>
then,<BR>
<BR>
I got the following output:<BR>
k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s)<BR>
k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s)<BR>
k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s)<BR>
k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<BR>
k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<BR>
k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<BR>
k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<BR>
k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<BR>
k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<BR>
k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s)<BR>
k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s)<BR>
k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s)<BR>
k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s)<BR>
k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s)<BR>
k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s)<BR>
k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s)<BR>
<BR>
---------------------------------------------------------------------<BR>
<BR>
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!!<BR>
<BR>
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 . <BR>
<BR>
Any thoughts ?<BR>
<BR>
<BR>
regards,<BR>
Jes.<BR>
<BR>
<BR>
On Thu, 12 Mar 2009 Berk Hess wrote :<BR>
><BR>
>Hi,<BR>
><BR>
>This is a very inefficient method for determining the viscosity.<BR>
>Also you need really perfect pressure fluctuations: NVT, shifted potentials,<BR>
>probably even double precision.<BR>
>There was a mail about this recently.<BR>
>There are better methods, have a look at:<BR>
>http://dx.doi.org/10.1063/1.1421362<BR>
><BR>
>Berk<BR>
><BR>
>Date: Thu, 12 Mar 2009 07:39:52 +0000<BR>
> From: jesbman@rediffmail.com<BR>
>To: gmx-users@gromacs.org<BR>
>Subject: Re: Re: [gmx-users] viscosity calculation using g_energy (3.3.3)<BR>
>CC:<BR>
><BR>
><BR>
>David,<BR>
><BR>
><BR>
><BR>
>Thanks for the quick reply.<BR>
><BR>
><BR>
><BR>
>Indeed I did as what you suggested- g_energy -f water.edr -vis test.xvg<BR>
><BR>
><BR>
><BR>
>The output file created includes three columns.<BR>
><BR>
><BR>
><BR>
>1. time ( ps) 2. shear viscosity (3) I assume it is bulk viscosity.<BR>
><BR>
><BR>
><BR>
>It seems, the unit given is cp. ( 1cp= 1* 10¯3 Pascal Second).<BR>
><BR>
><BR>
><BR>
>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.<BR>
><BR>
><BR>
><BR>
>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.<BR>
><BR>
><BR>
><BR>
>I post the first 20 lines of the output file.<BR>
><BR>
><BR>
><BR>
>-------------------------------------------------------------------<BR>
><BR>
><BR>
><BR>
># This file was created Thu Mar 12 16:20:09 2009<BR>
><BR>
># by the following command:<BR>
><BR>
># g_energy -f water.edr -vis test.xvg<BR>
><BR>
>#<BR>
><BR>
># g_energy is part of G R O M A C S:<BR>
><BR>
>#<BR>
><BR>
># GROup of MAchos and Cynical Suckers<BR>
><BR>
>#<BR>
><BR>
>@ title "Bulk Viscosity"<BR>
><BR>
>@ xaxis label "Time (ps)"<BR>
><BR>
>@ yaxis label "\8h\4 (cp)"<BR>
><BR>
>@TYPE xy<BR>
><BR>
>@ view 0.15, 0.15, 0.75, 0.85<BR>
><BR>
>@ legend on<BR>
><BR>
>@ legend box on<BR>
><BR>
>@ legend loctype view<BR>
><BR>
>@ legend 0.78, 0.8<BR>
><BR>
>@ legend length 2<BR>
><BR>
>@ s0 legend "Shear"<BR>
><BR>
>@ s1 legend "Bulk"<BR>
><BR>
> 1.99203 9.6633 96.3893<BR>
><BR>
> 3.98406 11.1625 98.1365<BR>
><BR>
> 5.9761 12.6631 99.838<BR>
><BR>
> 7.96813 13.4652 101.366<BR>
><BR>
> 9.96016 13.7012 100.249<BR>
><BR>
>-------------------------------------------------------------------------<BR>
</P>
<br><br>
<Table border=0 Width=644 Height=57 cellspacing=0 cellpadding=0 style='font-family:Verdana;font-size:11px;line-height:15px;'><TR><td><a href='http://adworks.rediff.com/cgi-bin/AdWorks/click.cgi/www.rediff.com/signature-home.htm/1050715198@Middle5/2705210_2676711/2680353/1?PARTNER=3&OAS_QUERY=null' target=new ><img src ='http://imadworks.rediff.com/cgi-bin/AdWorks/adimage.cgi/2705210_2676711/creative_2680353.gif' alt='Moneywiz' border=0></a></td></TR></Table>