GROMACS version: 2020.4

GROMACS modification: No

Hello all,

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.

Thanks,

Matthew