Finite-size corrections with Yeh-Hummer#

When performing molecular dynamics simulations with periodic boundary conditions, the calculated diffusion coefficients are affected by finite-size effects. The kinisi package includes tools to apply the Yeh-Hummer finite-size correction, which extrapolates the diffusion coefficient to the infinite system size limit.

In this tutorial, we will demonstrate how to use the YehHummer class to correct diffusion coefficients obtained from simulations with different box sizes.

Note

The Yeh-Hummer correction is based on the work by Yeh & Hummer (J. Phys. Chem. B 2004, 108, 15873-15879). The method assumes that the finite-size effects follow a linear relationship with the inverse box length.

[1]:

# This initialises the duecredit framework
%env DUECREDIT_ENABLE=yes
%env DUECREDIT_FILE=yehhummer.duecredit.p
import duecredit
duecredit.due.activate()

import matplotlib.pyplot as plt
import numpy as np
import scipp as sc

from kinisi.yeh_hummer import YehHummer

np.random.seed(42)

env: DUECREDIT_ENABLE=yes
env: DUECREDIT_FILE=yehhummer.duecredit.p

Setting up the data#

For this tutorial, we will use the TIP3P water diffusion data from the original Yeh & Hummer paper. This data represents diffusion coefficients calculated from molecular dynamics simulations of water at different box sizes.

For obtaining diffusion coefficients from your simulations, you can use the kinisi.analyze.DiffusionAnalyzer (from DiffusionAnalyzer.D) by following the steps outlined in the Diffusion coefficient from a VASP file tutorial.

[2]:

# TIP3P water data from Yeh & Hummer paper
box_lengths = np.array([18.58, 23.42, 29.51, 37.19, 46.86])  # Angstroms
D_values = np.array([4.884e-5, 5.123e-5, 5.315e-5, 5.466e-5, 5.590e-5])  # cm^2/s
D_errors = np.array([0.032e-5, 0.027e-5, 0.014e-5, 0.011e-5, 0.013e-5])  # cm^2/s

We need to organize this data into a scipp.DataArray with the appropriate coordinates and variances.

[3]:

td = sc.DataArray(
    data=sc.array(dims=['system'], values=D_values, variances=D_errors**2, unit='cm^2/s'),
    coords={'box_length': sc.Variable(dims=['system'], values=box_lengths, unit='angstrom')},
)
td

[3]:

scipp.DataArray (1009 Bytes)

Dimensions:
- system: 5
Coordinates: (1)
- box_length
  (system)
  float64
  Å
  18.58, 23.42, 29.51, 37.19, 46.86
```
Values:
array([18.58, 23.42, 29.51, 37.19, 46.86])
```

Data: (1)

(system)

float64

cm^2/s

4.884e-05, 5.123e-05, 5.315e-05, 5.466e-05, 5.590e-05

σ = 3.200e-07, 2.700e-07, 1.400e-07, 1.100e-07, 1.300e-07

Values:
array([4.884e-05, 5.123e-05, 5.315e-05, 5.466e-05, 5.590e-05])

Variances (σ²):
array([1.024e-13, 7.290e-14, 1.960e-14, 1.210e-14, 1.690e-14])

Creating the YehHummer object#

With our data prepared, we can now create a YehHummer object. The temperature is required for the viscosity calculation. You can optionally provide bounds for the fitting parameters.

[4]:

# Create YehHummer object with temperature 298 K
yh = YehHummer(td, temperature=sc.scalar(value=298, unit='K'))

/opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/scipy/optimize/_numdiff.py:710: RuntimeWarning: invalid value encountered in subtract
  df = [f_eval - f0 for f_eval in f_evals]

MCMC sampling for uncertainty estimation#

To obtain full probability distributions for the parameters, we can perform Markov chain Monte Carlo (MCMC) sampling:

[5]:

yh.mcmc(n_samples=500, n_walkers=16)

MCMC Sampling: 100%|██████████| 1000/1000 [00:03<00:00, 312.18it/s]

After MCMC sampling, we can visualize the probability distributions for the fitted parameters:

[6]:

from corner import corner

corner(
    np.array([i.values for i in yh.flatchain.values()]).T,
    labels=[' / '.join([k, str(v.unit)]) for k, v in yh.flatchain.items()],
)
plt.show()

Accessing the results#

The infinite-system diffusion coefficient and shear viscosity are now available as probability distributions:

[7]:

yh.D_infinite

[7]:

kinisi.Samples

(samples: 800)

float64

cm^2/s

(6.052+/-0.027)e-05

Values:
array([6.09479064e-05, 6.01949976e-05, 5.97237330e-05, 6.07667734e-05,
       6.01239400e-05, 6.06239664e-05, 6.02146726e-05, 6.02841872e-05,
       6.03109463e-05, 6.05577974e-05, 6.09685772e-05, 6.03494702e-05,
       6.02747140e-05, 5.98728887e-05, 6.06815820e-05, 5.99794007e-05,
       6.06730267e-05, 6.04925241e-05, 6.05791738e-05, 6.09483879e-05,
       6.01512642e-05, 6.04299763e-05, 6.05018043e-05, 6.04994090e-05,
       6.06302017e-05, 6.06951329e-05, 6.06086768e-05, 6.06740416e-05,
       6.04080745e-05, 6.02073283e-05, 6.10291826e-05, 6.04954332e-05,
       6.12051897e-05, 6.06181120e-05, 6.03170322e-05, 6.06713936e-05,
       6.02793278e-05, 6.08924994e-05, 6.03355921e-05, 6.03862513e-05,
       6.00414309e-05, 6.05922176e-05, 6.07530509e-05, 6.04110242e-05,
       6.07609832e-05, 6.01026114e-05, 6.10426136e-05, 6.03730503e-05,
       6.07833021e-05, 6.06768235e-05, 6.04454904e-05, 6.08213754e-05,
       6.11054254e-05, 6.08754613e-05, 6.05271763e-05, 6.04110814e-05,
       6.04574877e-05, 6.05893152e-05, 6.01904677e-05, 6.03270494e-05,
       6.05203402e-05, 6.04525154e-05, 6.06294428e-05, 6.04632470e-05,
       6.06208173e-05, 6.02729700e-05, 6.06424437e-05, 6.08788528e-05,
       6.06277208e-05, 6.10522623e-05, 6.07004347e-05, 6.03100225e-05,
       6.03718042e-05, 6.08606352e-05, 5.98679999e-05, 6.04244445e-05,
       6.02555598e-05, 5.99084559e-05, 6.08472003e-05, 6.06574993e-05,
       6.04760839e-05, 6.02428910e-05, 6.04006962e-05, 6.10537785e-05,
       6.03513672e-05, 6.04451817e-05, 6.08420799e-05, 6.00568294e-05,
       6.05666626e-05, 6.06793867e-05, 5.99018280e-05, 6.03392056e-05,
       6.01986675e-05, 6.04342640e-05, 6.06316650e-05, 6.07280369e-05,
       6.05172397e-05, 6.04481737e-05, 6.00745694e-05, 6.11556526e-05,
       6.05643851e-05, 5.99852822e-05, 6.07888002e-05, 6.00047309e-05,
       6.05113187e-05, 6.05213980e-05, 6.00860103e-05, 6.01061646e-05,
       6.08032839e-05, 6.07538063e-05, 6.04941997e-05, 6.05222862e-05,
       6.03750440e-05, 6.04374218e-05, 6.00152934e-05, 6.09304497e-05,
       6.04010149e-05, 6.03604715e-05, 6.07193569e-05, 6.00312771e-05,
       6.07024886e-05, 6.04267876e-05, 6.01240681e-05, 6.04756363e-05,
       6.07330315e-05, 6.07879433e-05, 6.03924624e-05, 6.01107384e-05,
       6.06177071e-05, 6.07129216e-05, 6.00674024e-05, 6.05704667e-05,
       6.05037407e-05, 6.04832262e-05, 6.09442300e-05, 6.05382272e-05,
       6.06854860e-05, 6.06057522e-05, 6.02225578e-05, 6.03637760e-05,
       6.02514309e-05, 6.11989069e-05, 6.02520582e-05, 6.01147554e-05,
       6.04721793e-05, 6.05005440e-05, 6.01918657e-05, 6.03314774e-05,
       6.02949101e-05, 6.06632087e-05, 6.10725268e-05, 6.05805043e-05,
       6.08432821e-05, 6.07027053e-05, 6.05530671e-05, 6.07512324e-05,
       6.01815723e-05, 6.09417684e-05, 6.05306156e-05, 5.98132878e-05,
       6.06881693e-05, 6.07268150e-05, 6.03003979e-05, 6.04766207e-05,
       6.06764892e-05, 6.06341626e-05, 6.07334981e-05, 5.97971566e-05,
       6.06948511e-05, 6.05842348e-05, 6.04228249e-05, 6.08449969e-05,
       6.05683223e-05, 6.08295698e-05, 6.05416040e-05, 5.98218696e-05,
       6.05557136e-05, 6.07447075e-05, 6.02500486e-05, 6.03642570e-05,
       6.02845281e-05, 6.08434115e-05, 6.06872809e-05, 6.02105325e-05,
       6.06171624e-05, 6.04926502e-05, 6.06365032e-05, 6.05756044e-05,
       6.06555018e-05, 6.12089748e-05, 6.06436294e-05, 6.06465643e-05,
       6.01701523e-05, 6.08685258e-05, 6.02552635e-05, 6.03604929e-05,
       6.07255086e-05, 6.06392496e-05, 6.08123117e-05, 5.97483527e-05,
       6.03576776e-05, 6.01969838e-05, 6.05048744e-05, 6.07410701e-05,
       6.07156922e-05, 6.07496259e-05, 6.06947822e-05, 6.09006127e-05,
       6.03638340e-05, 6.06618936e-05, 6.04636722e-05, 6.06788708e-05,
       6.10449108e-05, 6.05566593e-05, 6.09264295e-05, 6.02913825e-05,
       6.05525669e-05, 6.05127106e-05, 6.06363664e-05, 6.09349192e-05,
       6.07387171e-05, 6.04065866e-05, 6.05299423e-05, 6.06206574e-05,
       6.06070858e-05, 6.05364887e-05, 6.02189371e-05, 6.10960843e-05,
       6.09928406e-05, 6.02099523e-05, 6.05293309e-05, 6.02859894e-05,
       6.04690463e-05, 6.05387583e-05, 6.05342869e-05, 6.10880496e-05,
       6.06856690e-05, 6.01265863e-05, 6.03604084e-05, 6.07288464e-05,
       6.09164930e-05, 6.04885468e-05, 6.00573951e-05, 6.07550034e-05,
       6.08711864e-05, 6.06553665e-05, 6.07272981e-05, 6.03149930e-05,
       6.04980207e-05, 6.04472815e-05, 6.02391800e-05, 6.10669895e-05,
       6.10208555e-05, 6.05461970e-05, 6.04334612e-05, 6.08611037e-05,
       6.05282561e-05, 6.03094417e-05, 6.03971313e-05, 6.04749991e-05,
       6.08776309e-05, 6.05066970e-05, 6.06238824e-05, 6.02354693e-05,
       6.04902884e-05, 6.05701824e-05, 6.04751552e-05, 6.11215731e-05,
       6.08321138e-05, 6.05263934e-05, 6.04132592e-05, 6.06909474e-05,
       6.06210354e-05, 6.01363516e-05, 6.08640908e-05, 6.02745587e-05,
       6.05648748e-05, 6.03846671e-05, 6.07758954e-05, 6.02862558e-05,
       6.05560358e-05, 6.04425931e-05, 6.07873715e-05, 6.09928165e-05,
       6.05421959e-05, 6.01188353e-05, 6.03786038e-05, 6.06487853e-05,
       6.08570741e-05, 6.02509883e-05, 6.05674501e-05, 5.98752977e-05,
       6.03630038e-05, 6.05170826e-05, 6.05149162e-05, 6.02603835e-05,
       6.12285120e-05, 6.04934180e-05, 6.06149142e-05, 6.05576132e-05,
       6.03237902e-05, 6.00576788e-05, 6.01879232e-05, 6.02306636e-05,
       6.06964195e-05, 6.07180437e-05, 6.07732419e-05, 6.02069207e-05,
       6.06132271e-05, 6.09495483e-05, 6.03396914e-05, 6.01224980e-05,
       6.06017834e-05, 6.03240229e-05, 6.04833636e-05, 6.03243203e-05,
       6.01801394e-05, 6.03477893e-05, 6.05518511e-05, 6.07602989e-05,
       6.07833191e-05, 6.06752654e-05, 6.11695619e-05, 6.06274545e-05,
       6.03447279e-05, 6.07705506e-05, 6.02133957e-05, 6.02261222e-05,
       6.04503821e-05, 6.03143007e-05, 6.04953807e-05, 6.06236360e-05,
       6.03837646e-05, 6.07046223e-05, 6.05345656e-05, 6.06317071e-05,
       6.07839643e-05, 6.08770682e-05, 6.10334851e-05, 6.05374254e-05,
       6.07496721e-05, 6.05657810e-05, 6.01055330e-05, 6.01626546e-05,
       6.03690977e-05, 6.06678993e-05, 6.06472920e-05, 6.04211781e-05,
       6.05917930e-05, 6.02289735e-05, 6.07652177e-05, 6.05553247e-05,
       6.08117907e-05, 6.08368253e-05, 6.06302867e-05, 6.08123100e-05,
       6.05219948e-05, 6.06710814e-05, 6.00643837e-05, 6.02474590e-05,
       6.04615229e-05, 6.06811007e-05, 6.04172729e-05, 6.11063962e-05,
       6.03508663e-05, 6.09039231e-05, 6.06254891e-05, 6.04624287e-05,
       6.03897070e-05, 6.08498608e-05, 6.07280969e-05, 6.07905367e-05,
       6.07133856e-05, 6.04248001e-05, 6.05773988e-05, 6.03068857e-05,
       6.01821696e-05, 6.02868410e-05, 6.07016705e-05, 6.10247028e-05,
       6.02883288e-05, 6.06680243e-05, 6.04835157e-05, 6.08288487e-05,
       5.99084073e-05, 6.07372133e-05, 6.05877747e-05, 6.06019377e-05,
       6.04893171e-05, 6.06425119e-05, 6.01956890e-05, 6.07235513e-05,
       6.04250937e-05, 6.04921127e-05, 6.01084057e-05, 6.08874924e-05,
       6.03246796e-05, 6.06882212e-05, 6.03232233e-05, 6.05209208e-05,
       6.01322014e-05, 6.09455121e-05, 6.05588988e-05, 6.06752724e-05,
       6.02075706e-05, 6.12493492e-05, 6.01320445e-05, 6.05086059e-05,
       6.04431250e-05, 6.04822401e-05, 6.02027250e-05, 6.07652671e-05,
       6.06222673e-05, 6.04054671e-05, 6.05714573e-05, 6.06712678e-05,
       6.03796406e-05, 6.11034875e-05, 6.06083665e-05, 6.03854615e-05,
       6.03745017e-05, 6.10426096e-05, 6.02111295e-05, 6.01162518e-05,
       6.01186825e-05, 6.03195339e-05, 6.02450678e-05, 6.03786278e-05,
       6.06117530e-05, 6.03002742e-05, 6.09940838e-05, 6.07522821e-05,
       6.03156743e-05, 6.06435097e-05, 6.03183026e-05, 6.07228980e-05,
       6.04772548e-05, 6.06413746e-05, 5.99611087e-05, 6.02217417e-05,
       6.02380880e-05, 6.01414339e-05, 6.02336722e-05, 6.00516781e-05,
       6.03393965e-05, 6.02930872e-05, 6.04643980e-05, 6.05822690e-05,
       6.01787446e-05, 6.06719800e-05, 6.01857304e-05, 6.06317351e-05,
       6.07648938e-05, 6.03074813e-05, 6.02725328e-05, 6.06634520e-05,
       6.05731607e-05, 6.03642475e-05, 6.03452591e-05, 6.01737826e-05,
       6.05204796e-05, 6.07970605e-05, 6.06560738e-05, 6.07649466e-05,
       6.00464049e-05, 6.05096887e-05, 6.02862231e-05, 6.05754930e-05,
       6.07271198e-05, 6.08710384e-05, 5.98775902e-05, 6.03240234e-05,
       6.08615048e-05, 6.06728783e-05, 6.03462505e-05, 6.02451430e-05,
       6.10049504e-05, 6.00962171e-05, 6.04782601e-05, 6.05498880e-05,
       6.02374911e-05, 6.04238909e-05, 6.05318100e-05, 6.05909500e-05,
       6.07012577e-05, 6.05289731e-05, 6.03165703e-05, 6.04174805e-05,
       6.09594657e-05, 6.07864823e-05, 6.01797726e-05, 6.04812381e-05,
       6.05784413e-05, 6.06620187e-05, 6.05029855e-05, 6.09868737e-05,
       6.04777915e-05, 6.06083226e-05, 6.01524949e-05, 6.06981405e-05,
       6.06856398e-05, 6.02929256e-05, 6.03345706e-05, 6.06866696e-05,
       6.08308035e-05, 6.08847559e-05, 6.03285290e-05, 6.02604343e-05,
       6.05009865e-05, 6.07443000e-05, 6.06342825e-05, 6.07026116e-05,
       6.02792306e-05, 6.06626959e-05, 6.05432772e-05, 6.08327469e-05,
       6.07646350e-05, 6.01784825e-05, 6.06043237e-05, 6.06307584e-05,
       6.08294210e-05, 6.03624593e-05, 6.05060071e-05, 6.06561954e-05,
       6.05758664e-05, 6.02202090e-05, 6.05563806e-05, 6.11257732e-05,
       6.04877000e-05, 6.08921570e-05, 6.06482341e-05, 6.07003844e-05,
       6.08187661e-05, 6.00830270e-05, 6.03216504e-05, 6.06614666e-05,
       6.09440524e-05, 6.03248406e-05, 6.03946281e-05, 6.07057319e-05,
       6.04997087e-05, 6.02495253e-05, 6.04946451e-05, 6.03065604e-05,
       6.04972987e-05, 6.07662409e-05, 6.06956650e-05, 6.07826413e-05,
       6.08424736e-05, 6.08418380e-05, 6.06474716e-05, 6.07546221e-05,
       6.06720893e-05, 6.01024240e-05, 6.04031606e-05, 6.05076815e-05,
       6.01097603e-05, 6.01585880e-05, 6.03909945e-05, 6.03487178e-05,
       6.04819117e-05, 6.05697848e-05, 6.05994128e-05, 6.07216666e-05,
       6.04944573e-05, 6.11746819e-05, 6.05230280e-05, 6.06728913e-05,
       6.05048328e-05, 6.06473665e-05, 6.01522084e-05, 6.05706762e-05,
       6.00805239e-05, 6.05140870e-05, 6.04773812e-05, 6.04566310e-05,
       6.04816191e-05, 6.05597197e-05, 6.04737863e-05, 6.06279836e-05,
       6.04552505e-05, 6.07100799e-05, 6.06081803e-05, 6.07104367e-05,
       6.04184568e-05, 6.07759395e-05, 6.00651926e-05, 6.05894120e-05,
       6.06366303e-05, 6.07157129e-05, 6.03875916e-05, 6.03389133e-05,
       6.04977606e-05, 6.03923435e-05, 6.04749459e-05, 6.05688097e-05,
       6.08695105e-05, 6.06764633e-05, 6.04464977e-05, 6.04218459e-05,
       6.06034857e-05, 6.08028965e-05, 6.00612526e-05, 6.05611121e-05,
       6.03957919e-05, 6.11040370e-05, 6.03454685e-05, 6.04203178e-05,
       6.06135438e-05, 6.05289970e-05, 6.03358802e-05, 6.01757428e-05,
       6.06370941e-05, 6.08535502e-05, 6.05413002e-05, 6.05501599e-05,
       6.06980941e-05, 6.08534541e-05, 6.01132902e-05, 6.06095224e-05,
       6.03535116e-05, 6.06332824e-05, 6.03168763e-05, 6.02895456e-05,
       6.03884397e-05, 6.05399452e-05, 6.04793345e-05, 6.02698740e-05,
       6.04931821e-05, 6.07171743e-05, 6.03754944e-05, 6.09354286e-05,
       6.04092653e-05, 6.11033298e-05, 6.01088761e-05, 6.06222975e-05,
       6.06923722e-05, 6.06650919e-05, 6.01904965e-05, 6.02165044e-05,
       5.99970343e-05, 6.05268336e-05, 6.04700588e-05, 6.06679361e-05,
       6.04965507e-05, 6.04525396e-05, 6.03208034e-05, 6.05176740e-05,
       6.06632478e-05, 6.04399440e-05, 6.00807645e-05, 6.03884352e-05,
       6.07127206e-05, 6.03828836e-05, 6.04085136e-05, 6.01168034e-05,
       6.01556519e-05, 6.06304617e-05, 6.09269222e-05, 6.03732546e-05,
       6.01794835e-05, 6.02781185e-05, 6.03227345e-05, 6.06923449e-05,
       6.05155877e-05, 6.05371479e-05, 6.00174251e-05, 6.05456951e-05,
       6.03323284e-05, 6.04782321e-05, 6.02693848e-05, 6.01816949e-05,
       6.02740863e-05, 6.07084165e-05, 6.09509248e-05, 6.04683067e-05,
       6.03037251e-05, 6.00337975e-05, 6.02930649e-05, 6.06358379e-05,
       6.10729661e-05, 6.06301996e-05, 6.01854507e-05, 6.08619963e-05,
       6.04335322e-05, 6.03806457e-05, 6.02909171e-05, 6.01626669e-05,
       6.04836145e-05, 6.05939563e-05, 6.03251668e-05, 6.06630477e-05,
       6.04087756e-05, 6.02304984e-05, 6.05579352e-05, 6.06235705e-05,
       6.05630475e-05, 6.08700926e-05, 6.01735873e-05, 6.05336029e-05,
       6.01993803e-05, 6.04878640e-05, 6.05010306e-05, 6.04175528e-05,
       6.08585693e-05, 6.09342280e-05, 5.99564403e-05, 6.09584013e-05,
       6.09928557e-05, 6.02497007e-05, 6.05044633e-05, 6.06041636e-05,
       6.05913862e-05, 6.08073713e-05, 5.99877843e-05, 6.08012272e-05,
       5.98522146e-05, 6.03244295e-05, 6.05174238e-05, 6.06139183e-05,
       6.05074069e-05, 6.09807577e-05, 6.03201740e-05, 6.07421536e-05,
       6.12700129e-05, 6.03016641e-05, 6.06890118e-05, 6.03988908e-05,
       6.07138063e-05, 6.04230623e-05, 6.04090707e-05, 6.06001114e-05,
       6.01605536e-05, 6.03434106e-05, 6.05909121e-05, 6.04177712e-05,
       6.06219165e-05, 6.09244055e-05, 6.03751142e-05, 6.12309807e-05,
       6.10286519e-05, 6.02793629e-05, 6.06532418e-05, 6.02347475e-05,
       6.03807669e-05, 6.02058319e-05, 6.03206820e-05, 6.06981962e-05,
       6.05002236e-05, 6.02696227e-05, 6.05306683e-05, 6.04746153e-05,
       6.04446192e-05, 6.07821298e-05, 6.01733029e-05, 6.07450160e-05,
       6.11340408e-05, 6.04628230e-05, 6.07797031e-05, 6.01868021e-05,
       6.04445512e-05, 6.08545569e-05, 6.02884447e-05, 6.04663223e-05,
       6.06722376e-05, 6.04654753e-05, 6.02004407e-05, 6.06272681e-05,
       6.03989503e-05, 6.07398890e-05, 6.02445132e-05, 6.09940305e-05,
       6.12015327e-05, 6.03035326e-05, 6.11394340e-05, 6.00500638e-05,
       6.01374645e-05, 6.06345866e-05, 6.06870674e-05, 6.07538796e-05,
       6.06857027e-05, 6.02840696e-05, 6.04874811e-05, 6.04996190e-05,
       6.05005904e-05, 6.07257742e-05, 6.05660164e-05, 6.07990969e-05,
       6.08270449e-05, 6.08280588e-05, 6.07967109e-05, 6.03581613e-05,
       6.03901215e-05, 6.07770702e-05, 6.06840965e-05, 6.03170455e-05,
       6.08513113e-05, 6.07543177e-05, 6.02807627e-05, 6.03440002e-05])

[8]:

yh.shear_viscosity

[8]:

kinisi.Samples

(samples: 800)

float64

Pa*s

0.000285+/-0.000012

Values:
array([0.00027059, 0.00029675, 0.00032496, 0.00027656, 0.00030687,
       0.00028506, 0.00030379, 0.00029534, 0.00029323, 0.00028535,
       0.00027048, 0.00029329, 0.00029602, 0.00031436, 0.0002779 ,
       0.00031175, 0.00028015, 0.00028639, 0.00028366, 0.0002689 ,
       0.00030508, 0.00029259, 0.00028716, 0.00028135, 0.0002787 ,
       0.00027659, 0.00028124, 0.00027789, 0.00029059, 0.00030024,
       0.00026664, 0.0002846 , 0.00025604, 0.00028267, 0.00029206,
       0.00027979, 0.00029788, 0.00027007, 0.00029695, 0.0002874 ,
       0.00030602, 0.00027966, 0.00027565, 0.00028713, 0.00027677,
       0.00030573, 0.00026724, 0.00029338, 0.00027269, 0.00028337,
       0.0002863 , 0.00027143, 0.00026287, 0.00027336, 0.00028296,
       0.00028853, 0.00028474, 0.00028128, 0.00029781, 0.00029249,
       0.00028709, 0.00028961, 0.00028114, 0.00028829, 0.00028017,
       0.00029577, 0.00027979, 0.00027325, 0.0002767 , 0.0002656 ,
       0.00027979, 0.00029136, 0.0002906 , 0.00026553, 0.00030695,
       0.00028948, 0.00029956, 0.00031821, 0.00027352, 0.00028642,
       0.00028644, 0.00029317, 0.00029189, 0.0002623 , 0.00028945,
       0.00028921, 0.00027373, 0.00030474, 0.00028542, 0.00027079,
       0.00030616, 0.00029461, 0.00029964, 0.00028847, 0.00027562,
       0.00027904, 0.00028938, 0.00029072, 0.00030878, 0.00025724,
       0.00028154, 0.00030908, 0.00027596, 0.00031146, 0.00028827,
       0.00028215, 0.00030114, 0.00030284, 0.00026914, 0.00027316,
       0.00028146, 0.00028633, 0.00029504, 0.00029192, 0.00031249,
       0.00026954, 0.00028868, 0.00029433, 0.00027818, 0.00031136,
       0.00027574, 0.00029198, 0.00030249, 0.00028706, 0.00027187,
       0.00027135, 0.00028825, 0.0003073 , 0.00027833, 0.00027771,
       0.00030546, 0.00028278, 0.00028581, 0.00028317, 0.00026964,
       0.0002886 , 0.00027536, 0.00028198, 0.0002882 , 0.00029057,
       0.0002903 , 0.00025821, 0.00029245, 0.00030755, 0.00028006,
       0.0002834 , 0.00030116, 0.00029356, 0.00029886, 0.00027915,
       0.00026272, 0.00028452, 0.00027519, 0.00027774, 0.00027785,
       0.00027707, 0.00029287, 0.00026901, 0.00028204, 0.00032487,
       0.00027786, 0.00027242, 0.00029925, 0.00028452, 0.00027717,
       0.00027743, 0.00027489, 0.00031951, 0.00027989, 0.00028221,
       0.00028056, 0.00027304, 0.00028025, 0.00026969, 0.00028165,
       0.0003254 , 0.00028608, 0.00027645, 0.0002979 , 0.00028618,
       0.00029401, 0.0002726 , 0.000277  , 0.00029773, 0.00027874,
       0.00028485, 0.00027458, 0.0002834 , 0.00027752, 0.00025539,
       0.0002731 , 0.00028287, 0.00030408, 0.00027068, 0.00029807,
       0.00029717, 0.00027591, 0.00027869, 0.00027246, 0.00031906,
       0.00029208, 0.00029569, 0.00028648, 0.00027523, 0.00027564,
       0.0002727 , 0.00027185, 0.00027739, 0.00029542, 0.00027957,
       0.00028482, 0.00028246, 0.00026525, 0.00028097, 0.00026941,
       0.00029301, 0.00028262, 0.00028531, 0.00027978, 0.00026877,
       0.00026925, 0.00027984, 0.00028217, 0.00028169, 0.0002838 ,
       0.00028438, 0.00029616, 0.00026648, 0.0002676 , 0.00029686,
       0.00028229, 0.0002926 , 0.00028275, 0.00028576, 0.00028522,
       0.0002612 , 0.00027486, 0.00029758, 0.0002892 , 0.00028107,
       0.00026941, 0.00028576, 0.0003041 , 0.0002777 , 0.00027053,
       0.00027891, 0.00027467, 0.00029379, 0.00028355, 0.00028986,
       0.00029579, 0.00026175, 0.00026175, 0.00027962, 0.00028969,
       0.00027668, 0.00028502, 0.0002935 , 0.00029176, 0.00028671,
       0.00027327, 0.00028298, 0.00027679, 0.00029792, 0.000283  ,
       0.0002808 , 0.00028478, 0.00025946, 0.00027138, 0.0002792 ,
       0.00029224, 0.00028315, 0.00027847, 0.00030209, 0.00027166,
       0.00029571, 0.0002843 , 0.00028651, 0.00026986, 0.0002918 ,
       0.00028316, 0.00028392, 0.00027361, 0.00026459, 0.00028292,
       0.00030131, 0.00029333, 0.00027721, 0.00027319, 0.00029619,
       0.00028472, 0.00031257, 0.00029557, 0.00028311, 0.00028175,
       0.00029272, 0.00025892, 0.00027946, 0.00028079, 0.00028213,
       0.00028814, 0.0003031 , 0.00030227, 0.00029752, 0.00027784,
       0.00027811, 0.00027373, 0.00029825, 0.00028265, 0.00026902,
       0.00029142, 0.00029918, 0.00028022, 0.00028986, 0.00028426,
       0.0002876 , 0.00029515, 0.00029244, 0.0002866 , 0.00027773,
       0.00027534, 0.00027815, 0.00025894, 0.00027881, 0.00029116,
       0.00027419, 0.00030405, 0.00029546, 0.00028638, 0.00029247,
       0.00027666, 0.00027574, 0.00029092, 0.00028043, 0.0002887 ,
       0.00027969, 0.00027479, 0.00027466, 0.00026229, 0.00028336,
       0.00027422, 0.00028329, 0.00030539, 0.00029782, 0.00029094,
       0.00027888, 0.00027768, 0.00028299, 0.0002826 , 0.00029705,
       0.00027817, 0.0002856 , 0.00027127, 0.00027471, 0.0002778 ,
       0.00027755, 0.00028542, 0.00027912, 0.00030513, 0.00029243,
       0.00028528, 0.00027645, 0.00028842, 0.00026019, 0.00029612,
       0.00027048, 0.00028144, 0.00028618, 0.00029044, 0.0002693 ,
       0.00027441, 0.00027558, 0.00027766, 0.00028645, 0.00027727,
       0.0002933 , 0.0002975 , 0.00028964, 0.00027814, 0.00026135,
       0.0002972 , 0.00027654, 0.00028501, 0.00026845, 0.00031931,
       0.00027408, 0.00027931, 0.00028419, 0.000288  , 0.00027878,
       0.00029835, 0.00027555, 0.00028532, 0.00028149, 0.00030379,
       0.00027152, 0.00029552, 0.00027872, 0.00029243, 0.00028193,
       0.00030671, 0.00026611, 0.00028357, 0.00028287, 0.00030284,
       0.00025654, 0.00029959, 0.00028182, 0.00028957, 0.00028661,
       0.00030027, 0.00027583, 0.00028359, 0.00029176, 0.0002835 ,
       0.00027661, 0.00029691, 0.00026008, 0.00028625, 0.00029062,
       0.00029453, 0.00026519, 0.00029863, 0.00029608, 0.00029899,
       0.00029344, 0.00029886, 0.00028953, 0.00028233, 0.00029474,
       0.00026386, 0.00027796, 0.00029476, 0.00027898, 0.00029433,
       0.00027933, 0.00028572, 0.00028089, 0.00030837, 0.00029264,
       0.0002948 , 0.00029787, 0.00029767, 0.00030402, 0.00029398,
       0.0002925 , 0.0002861 , 0.00028514, 0.0002995 , 0.00027682,
       0.00030141, 0.0002794 , 0.00027496, 0.00029399, 0.00029393,
       0.0002753 , 0.00028127, 0.00028952, 0.00029469, 0.00029857,
       0.00028421, 0.00027417, 0.00027482, 0.00027727, 0.0003014 ,
       0.0002833 , 0.00029569, 0.00028285, 0.00027894, 0.0002727 ,
       0.0003166 , 0.00028887, 0.00027528, 0.00027671, 0.00029119,
       0.00029257, 0.00026592, 0.00030458, 0.00028131, 0.00028366,
       0.00029544, 0.00028473, 0.00028225, 0.00027862, 0.00027843,
       0.00028621, 0.00029486, 0.00028661, 0.00026874, 0.00027238,
       0.0002972 , 0.00028446, 0.00028266, 0.00027922, 0.00028435,
       0.00026695, 0.00028671, 0.00028034, 0.00029333, 0.00027325,
       0.00028006, 0.00029483, 0.00029496, 0.00027585, 0.00027523,
       0.00027091, 0.00029398, 0.00029176, 0.00028316, 0.00027437,
       0.00028098, 0.00027764, 0.00029661, 0.0002785 , 0.00027985,
       0.00027085, 0.00027187, 0.00029932, 0.00028271, 0.00027644,
       0.00027701, 0.00029068, 0.00028581, 0.00027561, 0.00028261,
       0.00029297, 0.00028371, 0.00026389, 0.00028908, 0.00027055,
       0.00027713, 0.00027623, 0.0002692 , 0.00030243, 0.00029499,
       0.00027368, 0.00026989, 0.00029318, 0.00029913, 0.00027366,
       0.00028414, 0.00029302, 0.00028203, 0.00029948, 0.00028516,
       0.00028121, 0.00027428, 0.00027387, 0.00027108, 0.00027315,
       0.00027969, 0.00026941, 0.00027672, 0.00029895, 0.0002901 ,
       0.00027827, 0.0003001 , 0.00029675, 0.00028392, 0.00029732,
       0.00028597, 0.00028753, 0.00027638, 0.00027787, 0.00028967,
       0.00025911, 0.0002825 , 0.00027584, 0.00028596, 0.00027876,
       0.00030021, 0.00028041, 0.00030196, 0.00028515, 0.0002839 ,
       0.00028713, 0.00028504, 0.00028468, 0.000278  , 0.00028079,
       0.00029082, 0.00027483, 0.00027783, 0.00027246, 0.00029275,
       0.0002722 , 0.00030444, 0.00028003, 0.00027854, 0.00027541,
       0.00028744, 0.00029842, 0.00028035, 0.00029346, 0.00028324,
       0.00028421, 0.00027241, 0.000276  , 0.00028927, 0.00028634,
       0.00027976, 0.00027122, 0.00029956, 0.00028161, 0.00028399,
       0.00026091, 0.00028992, 0.00028791, 0.00027883, 0.00028668,
       0.00028914, 0.00030501, 0.00028206, 0.00026699, 0.00028393,
       0.00028353, 0.00027634, 0.00027058, 0.00029749, 0.00028329,
       0.00028952, 0.00027668, 0.00029236, 0.00029288, 0.00028842,
       0.00028479, 0.00028614, 0.00029883, 0.00028914, 0.0002742 ,
       0.0002904 , 0.00026595, 0.00029034, 0.00025977, 0.00030122,
       0.00028193, 0.00027653, 0.0002744 , 0.00030331, 0.00029468,
       0.00030719, 0.00028317, 0.00028936, 0.00027974, 0.00029242,
       0.00028758, 0.00029303, 0.00028616, 0.00028254, 0.00028944,
       0.00030275, 0.0002938 , 0.00027525, 0.00028437, 0.00029264,
       0.00030387, 0.00030136, 0.0002751 , 0.00026844, 0.00029076,
       0.00030951, 0.0002945 , 0.00029329, 0.0002776 , 0.00028807,
       0.00029027, 0.00030784, 0.00028312, 0.00029576, 0.00028192,
       0.0002946 , 0.0003043 , 0.00029669, 0.00027491, 0.00026769,
       0.00028481, 0.00029499, 0.00030547, 0.0002899 , 0.00027857,
       0.0002646 , 0.00028239, 0.00030032, 0.00026911, 0.00029101,
       0.00028631, 0.00029216, 0.00030581, 0.00028831, 0.00028121,
       0.00029857, 0.00027909, 0.00028904, 0.00029337, 0.00028313,
       0.00028385, 0.0002851 , 0.00027426, 0.00030384, 0.00028538,
       0.00030039, 0.00028512, 0.00028347, 0.00028939, 0.00027084,
       0.00026464, 0.00031537, 0.00026901, 0.00026396, 0.00029159,
       0.00028907, 0.00028313, 0.0002827 , 0.00027454, 0.00031333,
       0.00027593, 0.00031684, 0.0002908 , 0.00028651, 0.00028014,
       0.00028448, 0.0002642 , 0.00029374, 0.00027626, 0.00025407,
       0.000292  , 0.00028117, 0.00028683, 0.0002784 , 0.00029136,
       0.00029134, 0.00028231, 0.00029967, 0.0002885 , 0.00027927,
       0.00028495, 0.00027864, 0.00026579, 0.00029368, 0.00025972,
       0.00026488, 0.00029347, 0.00028212, 0.00029321, 0.00029051,
       0.00029559, 0.00029371, 0.00027918, 0.00028268, 0.0002889 ,
       0.00028199, 0.00028398, 0.00028779, 0.00027244, 0.00030519,
       0.00027428, 0.00026033, 0.00028417, 0.00027606, 0.00029532,
       0.00029035, 0.00027135, 0.00029475, 0.00028907, 0.00027826,
       0.00028599, 0.00029713, 0.00027634, 0.00028989, 0.00027245,
       0.00029473, 0.00026407, 0.00025769, 0.00029326, 0.00026244,
       0.00029878, 0.00030229, 0.00028021, 0.00027983, 0.00027524,
       0.00027891, 0.00029702, 0.00028379, 0.00028083, 0.00029262,
       0.00027434, 0.00027923, 0.0002709 , 0.00027133, 0.00027461,
       0.00027226, 0.00028815, 0.00028965, 0.00027489, 0.00027502,
       0.00029244, 0.00027284, 0.00027458, 0.00029076, 0.00028911])

Plotting the results#

We can now plot the linear relationship between the diffusion coefficient and the inverse box length with credible intervals. For the estimation of the diffusion coefficient \(D_\infty\) we extrapolate to infinite box size \(\frac{1}{L} \to 0\):

[9]:

fig, ax = plt.subplots()

credible_intervals = [[16, 84], [2.5, 97.5], [0.15, 99.85]]
alpha = [0.6, 0.4, 0.2]

inv_L_data = 1 / yh.box_lengths.values
max_inv_L = np.max(inv_L_data)
inv_L_extended = np.linspace(0, max_inv_L * 1.1, 50)

ax.errorbar(
    inv_L_data,
    yh.diffusion.values,
    np.sqrt(yh.diffusion.variances),
    marker='o',
    ls='',
    color='k',
    zorder=10,
)

D_0_samples = yh.data_group['D_0'].values
slope_samples = yh.data_group['slope'].values
n_samples = len(D_0_samples)
predictions_extended = np.zeros((len(inv_L_extended), n_samples))

for i in range(n_samples):
    predictions_extended[:, i] = D_0_samples[i] - slope_samples[i] * inv_L_extended

for i, ci in enumerate(credible_intervals):
    ax.fill_between(
        inv_L_extended,
        *np.percentile(predictions_extended, ci, axis=1),
        alpha=alpha[i],
        color='#0173B2',
        lw=0,
    )

ax.set_xlabel(f'1 / L / {yh.box_lengths.unit}')
ax.set_ylabel(f'D  / {yh.diffusion.unit}')
plt.show()

Understanding the results#

The Yeh-Hummer correction provides two important quantities:

Infinite-system diffusion coefficient (:math:`D_{infty}`): This is the extrapolated diffusion coefficient for an infinite system.
Shear viscosity (:math:`eta`): This is estimated from the slope of the linear relationship:

\[D_{\text{PBC}} = D_\infty - \frac{k_B T \xi}{6\pi\eta L}\]

where:

\(D_{\text{PBC}}\) is the diffusion coefficient from the periodic simulation
\(D_\infty\) is the infinite-system diffusion coefficient
\(k_B\) is the Boltzmann constant
\(T\) is the temperature
\(\xi\) constant from Ewald summation (2.837297 for cubic boxes)
\(\eta\) is the shear viscosity
\(L\) is the box length

The duecredit framework can then be used to output the relevant references that should be cited.

[10]:

from duecredit import due

due.dump()


DueCredit Report:
- kinisi: Bayesian analysis of mass transport from molecular dynamics simulations / kinisi (v 2.1.0) [1]
  - Yeh-Hummer finite-size correction. / kinisi.yeh_hummer.YehHummer (v 2.1.0) [2]
- Scientific tools library / numpy (v 2.5.4) [3]

2 packages cited
1 module cited
0 functions cited

References
----------

[1] McCluskey, A.R. et al., 2024. kinisi: Bayesian analysis of mass transport from molecular dynamics simulations. Journal of Open Source Software, 9(94), p.5984.
[2] Yeh, I.-C. & Hummer, G., 2004. System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations with Periodic Boundary Conditions. The Journal of Physical Chemistry B, 108(40), pp.15873–15879.
[3] Van Der Walt, S., Colbert, S.C. & Varoquaux, G., 2011. The NumPy array: a structure for efficient numerical computation. Computing in Science & Engineering, 13(2), pp.22–30.