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: jesb...@rediffmail.com
>To: gmx-users@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
>
>-------------------------------------------------------------------------
_______________________________________________
gmx-users mailing list gmx-users@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-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/mailing_lists/users.php
