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: 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
_______________________________________________ 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