[gmx-users] viscosity calculation using g_energy (3.3.3)
JMandumpal
jesbman at rediffmail.com
Wed Mar 25 03:12:15 CET 2009
Dear Berk,
Thanks for the reply. I have read through your paper. I have still some doubts.
When used tcafś I got files called
tcaf_all.xvg tcaf_fit.xvg tcaf.xvg visc_k.xvg ( all default names).
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?
IF yes, what all other files ( tcf_all.xvg, tcaf_fit.xvg and tcaf.xvg) do?
expecting your reply,
Jestin
On Fri, 13 Mar 2009 Berk Hess wrote :
>
>Hi,
>
>I don't understand what you are actually doing now.
>You seem to be mixing multiple methods.
>
>First off all, I would use NPT for all methods, except the one that uses the pressure fluctuation.
>The pressure will have a large effect on the viscosity and if you run NVT you need to have
>exactly the right volume.
>
>If you use the cosine acceleration method, the 1/viscosity is printed in the energy file,
>g_energy will plot it for you.
>
>g_tcaf is only for use with an equilibrium simulation.
>If you read the paper, you will have seen an expression to extrapolate the k=0.
>
>Berk
>
>Date: Fri, 13 Mar 2009 07:33:30 +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:
>
>
> 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
>
> >
>
> >-------------------------------------------------------------------------
>
>
>
>
>
>
>
>
>_________________________________________________________________
>Express yourself instantly with MSN Messenger! Download today it's FREE!
>http://messenger.msn.click-url.com/go/onm00200471ave/direct/01/
>_______________________________________________
>gmx-users mailing list gmx-users at gromacs.org
>http://www.gromacs.org/mailman/listinfo/gmx-users
>Please search the archive at http://www.gromacs.org/search before posting!
>Please don't post (un)subscribe requests to the list. Use the
>www interface or send it to gmx-users-request at gromacs.org.
>Can't post? Read http://www.gromacs.org/mailing_lists/users.php
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090325/6129af19/attachment.html>
More information about the gromacs.org_gmx-users
mailing list