<P>
Dear Berk,<BR>
<BR>
Thanks for the reply. I have read through your paper. I have still some doubts.<BR>
<BR>
When used tcaf&#347; I got files called<BR>
<BR>
tcaf_all.xvg tcaf_fit.xvg tcaf.xvg visc_k.xvg ( all default names).<BR>
<BR>
I think, the file which I should use for fitting, using the formula eta(k) = eta (0) (1-ak²) + O(k^4) the viscoty ( viscosity - k vector plot) is visc_k.xvg. Am I right? <BR>
<BR>
IF yes, what all other files ( tcf_all.xvg, tcaf_fit.xvg and tcaf.xvg) do?<BR>
<BR>
expecting your reply,<BR>
Jestin<BR>
<BR>
On Fri, 13 Mar 2009 Berk Hess wrote :<BR>
><BR>
>Hi,<BR>
><BR>
>I don't understand what you are actually doing now.<BR>
>You seem to be mixing multiple methods.<BR>
><BR>
>First off all, I would use NPT for all methods, except the one that uses the pressure fluctuation.<BR>
>The pressure will have a large effect on the viscosity and if you run NVT you need to have<BR>
>exactly the right volume.<BR>
><BR>
>If you use the cosine acceleration method, the 1/viscosity is printed in the energy file,<BR>
>g_energy will plot it for you.<BR>
><BR>
>g_tcaf is only for use with an equilibrium simulation.<BR>
>If you read the paper, you will have seen an expression to extrapolate the k=0.<BR>
><BR>
>Berk<BR>
><BR>
>Date: Fri, 13 Mar 2009 07:33:30 +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>
> Dear Berk and David,<BR>
><BR>
><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>
><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>
><BR>
>Apart from that I added following options to my mdp file, where accelaration of 1A/ps² was given to the system.<BR>
><BR>
><BR>
><BR>
>;NON EQUILIBRIUM STUFF<BR>
><BR>
>acc_grps = system<BR>
><BR>
>accelerate = 0.1 0.0 0.0<BR>
><BR>
>cos_acceleration = 0.1<BR>
><BR>
><BR>
><BR>
>----------------------------------------------------------------<BR>
><BR>
><BR>
><BR>
>Moreover, I saved the trajectory in every 1ps ( so total 500 frames for a 500ps simulation)<BR>
><BR>
><BR>
><BR>
>then,<BR>
><BR>
><BR>
><BR>
>I got the following output:<BR>
><BR>
>k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s)<BR>
><BR>
>k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s)<BR>
><BR>
>k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s)<BR>
><BR>
>k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<BR>
><BR>
>k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<BR>
><BR>
>k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<BR>
><BR>
>k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<BR>
><BR>
>k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<BR>
><BR>
>k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<BR>
><BR>
>k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s)<BR>
><BR>
>k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s)<BR>
><BR>
>k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s)<BR>
><BR>
>k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s)<BR>
><BR>
>k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s)<BR>
><BR>
>k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s)<BR>
><BR>
>k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s)<BR>
><BR>
><BR>
><BR>
>---------------------------------------------------------------------<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>
><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>
><BR>
><BR>
>Any thoughts ?<BR>
><BR>
><BR>
><BR>
><BR>
><BR>
>regards,<BR>
><BR>
>Jes.<BR>
><BR>
><BR>
><BR>
><BR>
><BR>
>On Thu, 12 Mar 2009 Berk Hess wrote :<BR>
><BR>
> ><BR>
><BR>
> >Hi,<BR>
><BR>
> ><BR>
><BR>
> >This is a very inefficient method for determining the viscosity.<BR>
><BR>
> >Also you need really perfect pressure fluctuations: NVT, shifted potentials,<BR>
><BR>
> >probably even double precision.<BR>
><BR>
> >There was a mail about this recently.<BR>
><BR>
> >There are better methods, have a look at:<BR>
><BR>
> >http://dx.doi.org/10.1063/1.1421362<BR>
><BR>
> ><BR>
><BR>
> >Berk<BR>
><BR>
> ><BR>
><BR>
> >Date: Thu, 12 Mar 2009 07:39:52 +0000<BR>
><BR>
> > From: jesbman@rediffmail.com<BR>
><BR>
> >To: gmx-users@gromacs.org<BR>
><BR>
> >Subject: Re: Re: [gmx-users] viscosity calculation using g_energy (3.3.3)<BR>
><BR>
> >CC:<BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> >David,<BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> >Thanks for the quick reply.<BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> >Indeed I did as what you suggested- g_energy -f water.edr -vis test.xvg<BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> >The output file created includes three columns.<BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> >1. time ( ps) 2. shear viscosity (3) I assume it is bulk viscosity.<BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> >It seems, the unit given is cp. ( 1cp= 1* 10¯3 Pascal Second).<BR>
><BR>
> ><BR>
><BR>
> ><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>
> ><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>
> ><BR>
><BR>
> ><BR>
><BR>
> >I post the first 20 lines of the output file.<BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> >-------------------------------------------------------------------<BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> ><BR>
><BR>
> ># This file was created Thu Mar 12 16:20:09 2009<BR>
><BR>
> ><BR>
><BR>
> ># by the following command:<BR>
><BR>
> ><BR>
><BR>
> ># g_energy -f water.edr -vis test.xvg<BR>
><BR>
> ><BR>
><BR>
> >#<BR>
><BR>
> ><BR>
><BR>
> ># g_energy is part of G R O M A C S:<BR>
><BR>
> ><BR>
><BR>
> >#<BR>
><BR>
> ><BR>
><BR>
> ># GROup of MAchos and Cynical Suckers<BR>
><BR>
> ><BR>
><BR>
> >#<BR>
><BR>
> ><BR>
><BR>
> >@ title "Bulk Viscosity"<BR>
><BR>
> ><BR>
><BR>
> >@ xaxis label "Time (ps)"<BR>
><BR>
> ><BR>
><BR>
> >@ yaxis label "\8h\4 (cp)"<BR>
><BR>
> ><BR>
><BR>
> >@TYPE xy<BR>
><BR>
> ><BR>
><BR>
> >@ view 0.15, 0.15, 0.75, 0.85<BR>
><BR>
> ><BR>
><BR>
> >@ legend on<BR>
><BR>
> ><BR>
><BR>
> >@ legend box on<BR>
><BR>
> ><BR>
><BR>
> >@ legend loctype view<BR>
><BR>
> ><BR>
><BR>
> >@ legend 0.78, 0.8<BR>
><BR>
> ><BR>
><BR>
> >@ legend length 2<BR>
><BR>
> ><BR>
><BR>
> >@ s0 legend "Shear"<BR>
><BR>
> ><BR>
><BR>
> >@ s1 legend "Bulk"<BR>
><BR>
> ><BR>
><BR>
> > 1.99203 9.6633 96.3893<BR>
><BR>
> ><BR>
><BR>
> > 3.98406 11.1625 98.1365<BR>
><BR>
> ><BR>
><BR>
> > 5.9761 12.6631 99.838<BR>
><BR>
> ><BR>
><BR>
> > 7.96813 13.4652 101.366<BR>
><BR>
> ><BR>
><BR>
> > 9.96016 13.7012 100.249<BR>
><BR>
> ><BR>
><BR>
> >-------------------------------------------------------------------------<BR>
><BR>
><BR>
><BR>
><BR>
><BR>
><BR>
><BR>
><BR>
>_________________________________________________________________<BR>
>Express yourself instantly with MSN Messenger! Download today it's FREE!<BR>
>http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/<BR>
>_______________________________________________<BR>
>gmx-users mailing list gmx-users@gromacs.org<BR>
>http://www.gromacs.org/mailman/listinfo/gmx-users<BR>
>Please search the archive at http://www.gromacs.org/search before posting!<BR>
>Please don't post (un)subscribe requests to the list. Use the<BR>
>www interface or send it to gmx-users-request@gromacs.org.<BR>
>Can't post? Read http://www.gromacs.org/mailing_lists/users.php<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://sigads.rediff.com/RealMedia/ads/click_nx.ads/www.rediffmail.com/signatureline.htm@Middle?" target="_blank"><IMG SRC="http://sigads.rediff.com/RealMedia/ads/adstream_nx.ads/www.rediffmail.com/signatureline.htm@Middle"></A></td></TR></Table>