I am trying to determine the density, self-diffusion coefficient, and shear viscosity of a system of diethylene glycol (DEG) (HOCH2CH2OCH2CH2OH) w/ small amounts of water (e.g. 900 DEG molecules + 100 water molecules of either SPC/E or TIP4P/2005) at 328K. Protocol wise, I generate my system using gmx insert-molecules, energy minimize using steepest descent, run through 1 ns of NPT equilibration with v-rescale (tau-t = 1.0 ps) and Parrinello-Rahman (tau-p = 5.0 ps) (I skip NVT equilibration since the system doesn’t crash), and then run 100 ns of NVT production with v-rescale. Note that I find the density from 200-1000 ps of the NPT equilibration run (since density is converged by that point). For the self-diffusion coefficient and shear viscosity, I find the MSD over 5000 position frames and record energies every 20 fs during the NVT production run since the shear-stress ACF requires frequent recording of energies to result in an accurate value for the shear viscosity.
Density and self-diffusion cause me no problems since the MSD > (box length = 5 nm)^2 (see first picture below) and I am sampling from the diffusive regime. However, viscosities found via gmx energy -vis -evisco to provide Green-Kubo integral and Einstein-relation values for shear viscosity do not appear to converge as seen in the second and third graphs, respectively, even though the production simulation is run for 100 ns at a moderate temperature. (Note that the blue line in the Einstein-relation graph is the average of the 3 off-diagonal elements of the shear-stress tensor.)
Thus, my question is as follows: Why do I not achieve viscosity convergence with either the GK integral or Einstein-relation approach? As such, which portion of the graph do I take for analysis given that the value for shear viscosity fluctuates so much over the time-axis (which I am aware, is actually the integration time, not simulation time) and how do I decide this objectively? In addition, I don’t believe I am doing anything wrong, so if someone could point out what is being done wrong, I would greatly appreciate it.
Note that the experimental viscosity for this system is approximately 8.8 cP. In addition, I have read the GROMACS manual sections regarding viscosity and also have read Berk Hess’s 2002 paper “Determining the shear viscosity of model liquids from molecular dynamics simulations”. Thus, I am aware of the methods for simulating viscosity which include both EMD and NEMD techniques.
These EMD viscosity problems have caused me quite a headache for a long time now, so if there is any explanation to this lack or convergence or any fix to allow me to objectively determine the viscosity, please let me know. I would greatly appreciate it.
Note that I take a frame of average density from NPT equilibration as a gro file and use that as my input into NVT production to set the volume. I then regenerate the velocities at 328K, which I am aware, may not be the best for the system. However, to compensate, I ignore an initial 5% of the run.
If there is anything with the mdp file that may be the cause of these non-convergence behaviors, please let me know.
Then I assume you setup is correct. It is difficult to converge measurements that are based on fluctuations. Double your simulation time, maybe multiple times, and see if you can obtain a plateau for the first few ns.
All different methods should converge to the same answer when set up correctly. But the challenges are different. With NEMD you need a low enough shear rate to get the correct viscosity.
I’ve done much NEMD simulations and checked many papers to find the most suitable shear rate for viscosity calculation. As far as my simulation results are concerned, smaller the shear rate is, higher viscosity you would finally get. And in most cases, there are different ‘‘low enough shear rate’’ for different species when simulation viscosity match the experimental value, so it means i choose a best shear rate for every different substance to make the simulation result more accurate. And that seems unreasonable. EMD don’t give good results either.
I am still learning the technicalities of viscosity determination, so I guess I cannot provide much feedback, but I can say something for NEMD: “low enough shear rate” really means approaching zero shear rate; that is the way to make it consistent with equilibrium Linear-response-theory approaces.
So what I have done previously is to simulate for several values of shear, maybe even logarithmically spaced, and then extrapolate for shear rate = 0.
I have tried this way. I mean in some particular shear rate, simulation have a good accuracy. When extrapolate for shear rate equals to 0, viscosity get further away from the experimental value. So which result should I choose? The value closest to truth or the zero shear rate point? I am really confused for a long time.
I have been working mainly with water, so I am not sure about your case.
Are you sure the viscosity from MD is going to match the one from experiments anyway?
For water models in almost never the case, although it is also generally the other way around w.r.t. your case: MD under-predicts experimental values, sometimes by quite a lot.
So I believe you can compare NEMD to linear-response results or to other MD estimates in the literature; direct comparison with experiments could be a long shot.
Thank you so much.
I know it’s difficult to obtain the high-accuracy viscosity through MD now, but in my case it is the basis of my follow-up work, so i’m trying to find some way to make the result more acceptable. Now i want to valid if the different methods could converge to the same answer as hess said, and i wonder if there are any advises or experience you can give me on how long the simulation time should be set cause the EMD method’s convergence is slow.
The convergence time depends on the system. If I am not mistaken, DEG is 35 more viscous than water, so you need much longer simulations that for water.
the experimental viscosity is compared with shear viscoisty. so you can see for your system near equal to 10-12cp and x axis around 2000 ps viscosity value converges. so I think like that we can consider the value. You can check for other system also if this works. I am not sure about this logic make sense or not. Also want to listen from hess’s thought on this
