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: jesb...@rediffmail.com
To: gmx-users@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: 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

>

>-------------------------------------------------------------------------








_________________________________________________________________
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@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

Reply via email to