[gmx-users] Potential bug? Different results from different versions of gromacs with cos_acceleration
zephyrwuz at hotmail.com
Wed Apr 25 04:45:30 CEST 2012
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!
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the gromacs.org_gmx-users