[gmx-users] viscosity calculation using g_energy (3.3.3)
Berk Hess
gmx3 at hotmail.com
Fri Mar 13 11:59:32 CET 2009
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/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090313/294401fa/attachment.html>
More information about the gromacs.org_gmx-users
mailing list