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]:

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

from kinisi.yeh_hummer import YehHummer

np.random.seed(42)

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'))

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:06<00:00, 158.66it/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

Finite-size corrections with Yeh-Hummer#

Setting up the data#

Creating the YehHummer object#

MCMC sampling for uncertainty estimation#

Accessing the results#

Plotting the results#

Understanding the results#

This Page