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: [email protected]
To: [email protected]
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                                                                      
          
-- 
gmx-users mailing list    [email protected]
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at 
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to [email protected].
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

Reply via email to