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
I am currently performing high-throughput NEMD simulations to determine the shear viscosity of complex, heavy hydrocarbon mixtures (diesel fuels, predominantly C15-C25 alkanes and aromatics) using the periodic perturbation method (cos-acceleration). My goal is to extract the zero-shear viscosity to rank 500 different compositions.
I have carefully read Dr. Hess’s 2002 paper (J. Chem. Phys. 116, 209-217) regarding this method. While the method works perfectly for simple liquids like SPC water or LJ fluids at moderate shear rates, I am encountering strong non-Newtonian behavior due to the long relaxation times of my heavy molecules, and I would greatly appreciate your insights on my extrapolation strategy and the observed systematic errors.
1. System Setup & Initial Observations:
System: ~36,000 atoms (mixture of long-chain alkanes and polycyclic aromatics).
Force Field: CHARMM36 (all-atom, non-polarizable).
Thermostat: Changed from Nose-Hoover to v-rescale (\tau=0.5 ps, T=293.15 K). (I found that Nose-Hoover caused severe temperature resonance \pm 30 K with the macroscopic flow, which was perfectly resolved by v-rescale).
Box: Rectangular (L_z \approx 21 nm).
Experimental Viscosity:\approx 5.0 cP.
2. The Shear Thinning Issue & Extrapolation Strategy:
When applying an acceleration of \mathcal{A} = 0.025 \text{ nm/ps}^2, I obtained a highly converged but very low apparent viscosity of 2.40 cP. Realizing this is deep in the shear-thinning regime, I performed a multi-point gradient test to extrapolate to a \rightarrow 0:
Performing a linear extrapolation using the lowest three shear points (\mathcal{A} = 0.005, 0.010, 0.015) yields a zero-shear viscosity intercept of approximately 9.1 cP.
3. My Questions for the Community/Developers:
Is this extrapolation mathematically/physically sound within the GROMACS NEMD framework? Since I cannot directly simulate at “sufficiently low” shear rates (e.g., \mathcal{A} = 0.0001) due to the poor Signal-to-Noise Ratio (SNR) overwhelmed by thermal fluctuations, is fitting the apparent viscosities at moderate \mathcal{A} values and extrapolating to a=0 the standard and accepted approach for long-chain polymers/heavy oils?
Regarding the absolute error: The extrapolated zero-shear viscosity (~9.1 cP) is almost 1.8 times higher than the macroscopic experimental value (~5.0 cP). Is it safe to assume that this overestimation is a systematic limitation of the rigid, non-polarizable CHARMM36 force field (which tends to over-predict intermolecular friction/entanglement for dense heavy alkanes)? Interestingly, the apparent viscosity at \mathcal{A} = 0.010 (5.9 cP) matches the experiment well, but I suspect this is merely a coincidence where the force field overestimation cancels out the shear-thinning underestimation.
High-Throughput Feature Extraction: For ranking 500 different systems, would it be scientifically justifiable to just simulate all of them at \mathcal{A} = 0.025 (where the SNR is excellent and converges in 5 ns) and use these “high-shear apparent viscosities” as relative feature descriptors for a Machine Learning model, rather than attempting the expensive multi-point extrapolation for all 500 systems?
Thank you very much for your time and guidance. I look forward to your expert opinions.
The extrapolation procedure looks reasonable. The viscosity on the force field and dihedral potentials are likely most critical. Charmm would not be my first choice for simulating polymers. Maybe the parameters are good, but I would suggest to go through the literature for validations.
For the viscosity calculation I would suggest the simple Couette flow setup that now work properly in the 2026 (and 2025) release. This allows for much lower shear rates for the same system size.
Thank you so much for your invaluable feedback! It is a great relief to know that the extrapolation procedure is fundamentally sound.
Your comment regarding the limitations of CHARMM dihedral potentials for long polymers/alkanes was incredibly thought-provoking. To investigate whether the overestimation was purely a force-field issue or a methodological artifact, I applied the exact same NEMD cos-acceleration extrapolation procedure to a second diesel mixture (System 2, a slightly lighter composition).
The results were fascinating:
System 2 (Lighter): Extrapolated zero-shear viscosity = 4.13 cP vs. Experimental = 4.09 cP(Near perfect match!)
System 1 (Heavier): Extrapolated zero-shear viscosity = ~9.1 cP vs. Experimental = ~5.0 cP(Large overestimation).
Given this perfect control experiment, I have two quick follow-up questions:
1. Force Field Recommendation: Since the methodology proved highly accurate for the lighter system, it seems the error in System 1 is purely driven by the force field (likely the artificial stiffness/friction of long chains in CHARMM36, as you pointed out). If CHARMM36 would not be your first choice for such heavy alkane mixtures, what force field would you strongly recommend evaluating first? (e.g., OPLS-AA, L-OPLS, TraPPE, or GROMOS?)
2. Couette Flow in GROMACS 2022: Your suggestion to use the steady-state simple Couette flow is brilliant, as the uniform shear rate would drastically improve the SNR. However, our supercomputing cluster is currently running GROMACS 2022. You specifically mentioned that this setup now works properly in the 2025/2026 releases. Does this imply that the Couette flow implementation (e.g., Lees-Edwards boundary conditions combined with thermostatting) has known bugs or stability issues in the 2022 release? Should I strictly avoid using deform for viscosity calculations in version 2022?
Thank you again for your time and continuous dedication to the community!
