Possible error in transport analysis?

[orionarcher]

I am working on a benchmark of viscosity. As part of this, I am comparing viscosity values of common electrolyte solvents to references in the literature. Strangely, I am getting an inverse trend, shown below. The values are in mPa s, the only data transformation I did was to change the units from Pa s to mPa s. Reference values are on the y axis.

As an experiment, I tried comparing the calculated data against the inverse of the reference data and the trend looks much better.

Obviously, there is a lot that could be causing this on the MD side but the workflow is reasonably rigorous. The viscosity is calculated from a 2 ns production run following 1.5 ns NPT and annealing, each data point is averaged over 5 replicates, linear fit window from 0.4 ns to 1.6 ns, 9 nm nonbonded cutoff. The fact that the order of magnitude is correct and the trend is wrong is very surprising, the opposite of what I'd expect. @hmacdope, I'm curious if you have any thoughts here? I've scoured the code and can't find any mistakes.

The non-log-log plots:

[hmacdope]

@orionarcher curious, I think your instincts are right that this is most likely a systematic error of some kind.

1. Did you use unwrapped trajectories? Worth trying with unwrapped positions.

2. Only thing I could spot from reading the code is that the reference we used for the implementation: https://iopscience.iop.org/article/10.1088/1742-6596/653/1/012106/pdf (eq.6 ) disagrees with https://livecomsjournal.org/index.php/livecoms/article/view/v1i1e6324/937 (Table 2), with the position of V in the prefactor (although included in subtext). I could be missing something, but https://arxiv.org/pdf/cond-mat/0701253 also agrees with Kirova & Norman, perhaps with $V_{inf}$ we need to remove. Unsure.

3. It may pay to do an exact reproduction of an LJ fluid, e.g from Kirova & Norman and see if we get the same numerical result. This should be pretty easy in LAMMPS.

4. We seem to be missing a velocity element and also not integrating correctly over the flux. See Table 2 Maginn.

5. We should check against GK estimate from stress tensor.

[hmacdope]

If you are OK sharing the data also, let me know and I can take a look.

[orionarcher]

0. I can share the data. It's taking more time than expected to get it into a suitable format but should be done tomorrow. My analysis scripts are annoyingly bespoke. The test data that I added to this repo is from the same dataset, if you'd like something to toy around with in the mean time.

1. I did use unwrapped trajectories. I thought that would be the problem but it doesn't appear to be.

2, 3, 4. I won't have a change to dig into the numerics until this weekend, at first glance I agree something looks off. Nothing would clearly generate the inverse relation that is appearing though. Very odd.

5. I reported every 2000 fs which, according to Maginn, is not enough to construct an accurate velocity autocorrelation. I'll try to run a trajectory over the weekend that includes more frequent reporting of velocities.

[rohithsrinivaas]

1. I did use unwrapped trajectories. I thought that would be the problem but it doesn't appear to be.

For Viscosity using the Helfand method, the wrapped coordinates needs to be used. Additionally two more terms are needed for the Helfand moment. These two terms are for the particles near the boundary. The details of the two terms are provided in https://arxiv.org/pdf/cond-mat/0510445 in pg. 53 eq. (2.67)

2, 3, 4. I won't have a change to dig into the numerics until this weekend, at first glance I agree something looks off. Nothing would clearly generate the inverse relation that is appearing though. Very odd.

The current code actually computes based on the Mc Quarrie's formula which is incorrect.

transport-analysis/transport_analysis/viscosity.py

Line 222 in 1199bfa

sq_diff = np.square(diff).mean(axis=-1)

The text from https://arxiv.org/pdf/cond-mat/0701253 clearly mentions this.