Viscosity calculation output?

GROMACS version: 2024.2
GROMACS modification: No

Hello ! I am trying to use the gmx energy [...] -vis -evisco Einstein viscosity analysis to calculate the viscosity of pure liquid succinonitrile, but I do not understand the units or how to use the resulting evisco.xvg. I have two questions, and then am attaching my mdp and plot of evisco.xvg for reference.

I ran 300 ns with 2000 SN molecules (20,000 atoms) using the OPLS forcefield. The viscosity of SN should be about 2 cP or 0.002 kg/m-sec according to experiments. Also, I have tried using nstenergy 1000 and also nstenergy 5, but got the same results, just with less noise.

Firstly, the evisco.xvg says the title is “shear viscosity using Einstein relation”, but what are the units in the “yaxis”: (kg m\S-1\N s\S-1\N)? The only capital S unit I know is Siemen, but there shouldn’t be any conductivity/charge measurements from my understanding… Are these columns truly the final viscosity prediction up to that time? or some intermediate term to be used in the Einstein relation?

Second, I noticed there are 5 columns in the evisco.xvg file, the first is clearly time, whereas the last is the average of the middle 3. Is this the viscosity calculated using integrals up to that time? or is this only the off-diagonal pressure terms? I ask because my average value of the xvg continues to grow linearly with time, even after 50 ns (post window averaging, I assume, in the xvg).

Should I expect the values of this xvg to plateau at longer simulation timescales?

Please let me know if you have any suggestions or require additional information to help. Thanks!

Pablo

SN_viscosity_evisco

nvt_prod_visc_einstein2.mdp (781 Bytes)

I’ll just answer your first question as I’m not an expert at the viscosity calculation at all. The yaxis label is written according to xmgrace format. If you open the xvg file in, e.g., xmgrace it should be rendered properly. I expect it to be kg m^-1 s^-1 with the ^-1 in superscript.

Hi,

The 4 columns after time corresponds to the three components of the off-diagonal stresses, + their average. This should be documented better.

Which file(s) are you looking at? I would expect the data in eviscoi.xvg to grow more or less linearly, while the output of evisco.xvg to stabilize (with possibly large fluctuations). The plot you’re showing seems a mix of the two.

Thank you both for the prompt replies!

Michele, your behavior is also what I would expect. Unfortunately, my evisco.xvg output is linear with time (except for one of the three off-diagonal stresses (blue)), while the output of the eviscoi.xvg is quadratic with time (since you take the linear evisco.xvg and integrate with respect to time, except for that first, blue, off-diagonal stress which was constant and is now linear).

Actually, if the eviscoi.xvg should increase linearly with time, but you have to square the integral before taking the slope (according to the gromacs documentation (Viscosity calculation - GROMACS 2024.3 documentation) equation 410), how would you find the slope? Because then the square would not be linear with time?

Regardless, I will look at using different force field parameters because I won’t be able to calculate a linear slope of my current eviscoi.xvg curves since they are not linear with time. Do you see anything else in the .mdb file I attached to my previous comment you would suggest changing when running a viscosity simulation? This is my first time calculating this property.

Thanks again,
Pablo

SN_viscosity_eviscoi

Are you deforming the box or applying an external force?
The only way I can explain this behaviour is that 2 out of 3 off-diagonal stress components are not zero on average (as a result of forcing/deformations or the system being very out of equilibrium).

Wow… after checking my equilibration .mdp files, I realized I had a typo and was doing my final stage of NPT equilibration before switching to NVT for production at 100 bars (instead of 1 bar). Therefore the system is most likely “very out of equilibrium”. I will rerun my equilibration and production runs and report back – in about a week.

Thank you for the time to offer suggestions!