<html>
<head>
<style>
.hmmessage P
{
margin:0px;
padding:0px
}
body.hmmessage
{
font-size: 10pt;
font-family:Verdana
}
</style>
</head>
<body class='hmmessage'>
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><hr id="stopSpelling">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>
Dear Berk and David,<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>
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>
Apart from that I added following options to my mdp file, where accelaration of 1A/psē was given to the system.<br>
<br>
;NON EQUILIBRIUM STUFF<br>
acc_grps = system<br>
accelerate = 0.1 0.0 0.0<br>
cos_acceleration = 0.1<br>
<br>
----------------------------------------------------------------<br>
<br>
Moreover, I saved the trajectory in every 1ps ( so total 500 frames for a 500ps simulation)<br>
<br>
then,<br>
<br>
I got the following output:<br>
k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s)<br>
k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s)<br>
k 1.593 tau 1.000 eta 0.09835 10^-3 kg/(m s)<br>
k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<br>
k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<br>
k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<br>
k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<br>
k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<br>
k 2.252 tau 1.000 eta 0.04917 10^-3 kg/(m s)<br>
k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s)<br>
k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s)<br>
k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s)<br>
k 2.759 tau 1.000 eta 0.03278 10^-3 kg/(m s)<br>
k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s)<br>
k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s)<br>
k 3.185 tau 1.000 eta 0.02459 10^-3 kg/(m s)<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>
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>
Any thoughts ?<br>
<br>
<br>
regards,<br>
Jes.<br>
<br>
<br>
On Thu, 12 Mar 2009 Berk Hess wrote :<br>
><br>
>Hi,<br>
><br>
>This is a very inefficient method for determining the viscosity.<br>
>Also you need really perfect pressure fluctuations: NVT, shifted potentials,<br>
>probably even double precision.<br>
>There was a mail about this recently.<br>
>There are better methods, have a look at:<br>
>http://dx.doi.org/10.1063/1.1421362<br>
><br>
>Berk<br>
><br>
>Date: Thu, 12 Mar 2009 07:39:52 +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>
>David,<br>
><br>
><br>
><br>
>Thanks for the quick reply.<br>
><br>
><br>
><br>
>Indeed I did as what you suggested- g_energy -f water.edr -vis test.xvg<br>
><br>
><br>
><br>
>The output file created includes three columns.<br>
><br>
><br>
><br>
>1. time ( ps) 2. shear viscosity (3) I assume it is bulk viscosity.<br>
><br>
><br>
><br>
>It seems, the unit given is cp. ( 1cp= 1* 10Ŋ3 Pascal Second).<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>
>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>
>I post the first 20 lines of the output file.<br>
><br>
><br>
><br>
>-------------------------------------------------------------------<br>
><br>
><br>
><br>
># This file was created Thu Mar 12 16:20:09 2009<br>
><br>
># by the following command:<br>
><br>
># g_energy -f water.edr -vis test.xvg<br>
><br>
>#<br>
><br>
># g_energy is part of G R O M A C S:<br>
><br>
>#<br>
><br>
># GROup of MAchos and Cynical Suckers<br>
><br>
>#<br>
><br>
>@ title "Bulk Viscosity"<br>
><br>
>@ xaxis label "Time (ps)"<br>
><br>
>@ yaxis label "\8h\4 (cp)"<br>
><br>
>@TYPE xy<br>
><br>
>@ view 0.15, 0.15, 0.75, 0.85<br>
><br>
>@ legend on<br>
><br>
>@ legend box on<br>
><br>
>@ legend loctype view<br>
><br>
>@ legend 0.78, 0.8<br>
><br>
>@ legend length 2<br>
><br>
>@ s0 legend "Shear"<br>
><br>
>@ s1 legend "Bulk"<br>
><br>
> 1.99203 9.6633 96.3893<br>
><br>
> 3.98406 11.1625 98.1365<br>
><br>
> 5.9761 12.6631 99.838<br>
><br>
> 7.96813 13.4652 101.366<br>
><br>
> 9.96016 13.7012 100.249<br>
><br>
>-------------------------------------------------------------------------<br>
<BR>
<br><br>
<table style="font-family: Verdana; font-size: 11px; line-height: 15px;" border="0" cellpadding="0" cellspacing="0" height="57" width="644"><tbody><tr><td><a href="http://adworks.rediff.com/cgi-bin/AdWorks/click.cgi/www.rediff.com/signature-home.htm/1050715198@Middle5/2705210_2676711/2680353/1?PARTNER=3&OAS_QUERY=null"><img src="http://imadworks.rediff.com/cgi-bin/AdWorks/adimage.cgi/2705210_2676711/creative_2680353.gif" alt="Moneywiz" border="0"></a></td></tr></tbody></table><br /><hr />Express yourself instantly with MSN Messenger! <a href='http://clk.atdmt.com/AVE/go/onm00200471ave/direct/01/' target='_new'>MSN Messenger</a></body>
</html>