[gmx-users] RE: Potential bug? Different results from different versions of gromacs with cos_acceleration

[Zhe Wu]

Just a follow up on this problem.
I think it is known that there is some bug in the 1/viscosity calculation in Gromacs version of 4.5.3. As shown in the following posts in the mailing list:1. lists.gromacs.org/pepermail/gmx-users/2011-January/057596.html2. lists.gromacs.org/pepermail/gmx-users/2011-March/059475.html
But both of them report the g_energy gives 1/viscosity as zero in the outcome, which is different from the case I reported here. Now with Gromacs 4.5.5, the previous bug could have been fixed. As stated in the release notes for 4.5.4, it says "Fixed incorrect cosine viscosity output". However, the output now is not zero but a number (viscosity) far too high compared to results from a different calculation method and an older version (3.3.3). Could it be a new bug in Gromacs 4.5.5?
One more details in the calculation with the periodic perturbation method: Since I am running NEMD and the method with Green-Kubo relation is for equilibrium MD, I just use g_energy without "-vis".
Thanks,Zhe

From: zephyrwuz at hotmail.com
To: gmx-users at gromacs.org
Subject: Potential bug? Different results from different versions of gromacs with cos_acceleration
Date: Wed, 25 Apr 2012 10:45:30 +0800





Dear Gromacs users and developers,
I found Gromacs 4.5.5 problematic in NEMD calculations for viscosity. As a simple benchmark, I did some calculations for water (TIP5P) and the results are confusing.
I applied two methods for the calculation in 300K, 512 waters: 1. momentum fluctuations, 2. periodic pertubation method, with Gromacs 4.5.5 double precision, each for 500 ps. (I understand it's short but it is just for benchmark and a longer calculation won't change the result qualitatively.) The first one is an equilibrium simulation (saving coordinates and velocity every 10 fs) and gives results as eta = 0.564 X 10^(-3) kg/(m * s), from the fitting of eta(k)=eta*(1-A*k^2) from g_tcaf, which is reasonable compared to experimental results and TIP5P results from earlier work. However, with periodic pertubation method, which is NEMD with cos_acceleration = 0.025 nm*ps^(-2) according to the original paper for SPC water, the result from g_energy gives 1/Viscosity as 103.276 m*s/kg. That means the viscosity is 9.68 X 10^(-3) kg/(m * s), which is an order of magnitude higher compared to the previous result! (Such result is reproduced for larger system, different cos_acceleration ( 
 < 1 nm*ps^(-2)) and longer runs.)
As a next step, I went into the code in mdlib and found that the code is relatively old. (I searched the mailing list and found most posts about such method is around 2007, and the ones around 2010 mentioning the similar observation but not responded.) So I used Gromacs 3.3.3 with the same inputs as previous for NEMD and performed the simulation again. Now with g_energy, 1/Viscosity is 2080.94 m*s/kg. (It doesn't mention specifically in the output for the unit, as stated as SI unit. So it should be m*s/kg, right?) Then the viscosity is 0.481 X 10^(-3) kg/(m * s), and qualitatively agrees with the previous results.
The original paper: http://jcp.aip.org/resource/1/jcpsa6/v116/i1/p209_s1
So the question is: Am I doing anything wrong? Or is the unit for viscosity in 3.3.3 different from 4.5.5? Or is NEMD just not a good method to calculate viscosity? (For the last question, it seems to be an open field for research.) I have to admit that the fluctuation of the results is very large (may due to the short simulation with small box) but I don't think it will change the results qualitatively.
Any suggestions will be really appreciated. Thank you in advanced for your help!
Thanks,Zhe 		 	   		   		 	   		  
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20120425/ee222784/attachment.html>