Kirkwood-Buff Integrator

Computes Kirkwood-Buff integrals (KBIs) from radial distribution functions (RDFs).

The KBI calculation typically suffers from finite-size effects because the RDF converges slowly toward the bulk value and is constrained by the simulation box. This class provides several correction schemes to recover the thermodynamic limit.

Default Behavior (weight_type='geometric')

By default, the rigorous framework outlined by Simon (2022) is employed to address finite-size effects:

  • RDF Correction: The RDF is adjusted for molecule excess/depletion (Ganguly and van der Vegt (2013)), ensuring the integral correctly reflects bulk density behavior.

  • Finite-Volume Correction: Corrects finite-volume effects by applying an analytical hypersphere model, compensating for volume loss at large radial distances (Krüger et al. (2013)).

  • Extrapolation: The integral is extrapolated to the thermodynamic limit, accounting for box-size dependence (Dawass, Krüger, et al. (2020)).

Correction Methods (weight_type parameter)

The weight_type parameter selects the specific weighting function \(w(r)\) or \(u_k(r)\) applied during integration. The 'geometric' option triggers the default behavior described above.

  • 'none': Uses the uncorrected RDF. Suitable for validation or systems where finite-size effects are negligible.

  • 'u0': Simple truncation of the integral at \(r=L\). Generally converges poorly as it lacks finite-volume corrections.

  • 'u1': First-order Taylor expansion in 1/L of the exact hypersphere weight function, \(w(r)\). Provides a more accurate estimate of \(G^\infty\) than the uncorrected u0 approach.

  • 'u2': Uses the exact hypersphere weight function with an analytical approximation of the surface term. This provides the most robust estimate of \(G^\infty\) among the non-extrapolated methods.

  • 'geometric': The most rigorous method. Uses the geometric weight function \(w(r)\) and performs explicit linear extrapolation to the thermodynamic limit.

Note

Linear extrapolation to the thermodynamic limit is performed only for weight_type='geometric'. The other weight functions (u0-u2) serve as direct estimates of \(G^\infty\) as demonstrated by Krüger and Vlugt (2018).

class kbkit.kbi.integrator.KBIntegrator(r: ndarray, gr: ndarray, n_ref: int, box_volume: float, delta: int, weight_type: Literal['none', 'u0', 'u1', 'u2', 'geometric'] = 'geometric', errors: Literal['raise', 'warn', 'ignore'] = 'raise', force: bool = False)[source]

Bases: object

Class to compute the Kirkwood-Buff Integrals (KBI) from RDF data.

Parameters:

  • r (np.ndarray) – Radial distance array (nm).

  • gr (np.ndarray) – Radial distribution function values.

  • box_volume (float) – Averaged simulation box volume (nm³).

  • delta (int) – Kronecker delta for RDF molecules.

  • weight_type (str, optional) – Type of weight function for finite-volume corrections. Options: (‘none’,’u0’,’u1’,’u2’,’geometric’)

  • errors (Literal["raise", "warn", "ignore"], optional) – Error mode for handling KBIConvergenceErrors. Only for weight_type='geometric'.

  • force (bool, optional) – If KBIConvergenceError is raised, warns user and returns KBI for weight_type='u2'. Only for weight_type='geometric'.

classmethod from_rdf(rdf: str | Path | RdfParser, system_properties: SystemProperties, weight_type: Literal['none', 'u0', 'u1', 'u2', 'geometric'] = 'geometric', errors: Literal['raise', 'warn', 'ignore'] = 'raise', force: bool = False) → KBIntegrator[source]

Create KBI object from RdfParser and SystemProperties.

Automatically extracts volume and molecule_count from system_properties.

Parameters:

  • rdf (str | Path | RdfParser) – RdfParser object or file containing rdf data. Supported filetypes: (‘.xvg’,’.txt’,’.csv’,’.xlsx’)

  • system_properties (SystemProperties) – System properties containing volume and topology.

  • weight_type (str, optional) – Type of weight function for finite-volume corrections. Options: (‘none’,’u0’,’u1’,’u2’,’geometric’)

  • errors (Literal["raise", "warn", "ignore"], optional) – Error mode for handling KBIConvergenceErrors. Only for weight_type='geometric'.

  • force (bool, optional) – If KBIConvergenceError is raised, warns user and returns KBI for weight_type='u2'. Only for weight_type='geometric'.

Returns:

Initialized integrator.

Return type:

KBI

property gr_vdv: ndarray

Compute the corrected radial distribution function, accounting for excess/depletion of molecule \(i\) around molecule \(j\) due to a finite number of molecules, based on the approach by Ganguly and Van der Vegt (2013).

Returns:

Factor for correcting RDF so density follows the bulk density of the system.

Return type:

np.ndarray

Notes

The correction is calculated as,

\[g^{vdV}(r) = g(r) \frac{N_j f(r)}{N_j f(r) - \Delta N_j(r) - \delta_{ij}}\]
\[f(r) = 1 - \frac{\frac{4}{3} \pi r^3}{\langle V \rangle}\]
\[\rho_j = \frac{N_j}{\langle V \rangle}\]
\[\Delta N_j(r) = \rho_j \int_0^{r_{max}} 4 \pi r^2 \bigl(g(r) - 1 \bigr) dr\]
where:

  • \(r\) is the distance

  • \(\langle V \rangle\) is the box volume

  • \(N_j\) is the number of particles of type ( j )

  • \(g(r)\) is the radial distribution function directly from simulation

  • \(\delta_{ij}\) is a kronecker delta

Note

The cumulative integral \(\Delta N_j(r)\) is approximated numerically using the trapezoidal rule.

compute_hr(weight_type: str) → ndarray[source]

Total correlation function, \(h(r) = g(r) - 1\).

If weight_type='none',

\[h(r) = g^{NpT}(r) - 1,\]

where \(g^{NpT}(r)\) is the uncorrected radial distribution function.

For any other weight_type specified, the Ganguly and Van der Vegt (2013) correction is applied,

\[h(r) = g^{vdV}(r) - 1.\]
property hr: ndarray

Total correlation function, \(h(r) = g(r) - 1\) based on the weight_type intitalized in class.

Type:

np.ndarray

geometric_weight(r: ndarray) → ndarray[source]

Geometric weighting function for 3D hypersphere.

Correct KBI for finite volumes with an exact analytically correct form for hyperspheres, based on the method described by Krüger et al. (2013).

Parameters:

r (np.ndarray) – Radial distance for weight function.

Returns:

Geometric weight function for hyperspheres.

Return type:

np.ndarray

Notes

The geometric weight function, \(w(r)\), is defined as:

\[w(r) = 4 \pi r^2 \left( 1 - \frac{3}{2}x + \frac{1}{2}x^3 \right),\]

where \(x\) is the dimensionless distance, \(x = r/L\).

u0_weight(r: ndarray) → ndarray[source]

Simple truncation of \(r = L\).

This weighting function is well-known to converge badly, it lacks any finite-volume corrections [Krüger and Vlugt (2018)].

Parameters:

r (np.ndarray) – Radial distance for weight function.

Returns:

Weight function.

Return type:

np.ndarray

Notes

The weight function, \(u_0(r)\), is defined as:

\[u_0(r) = 4 \pi r^2\]
where:

  • \(r\) is the radial distance

u1_weight(r: ndarray) → ndarray[source]

First order Taylor series expansion in \(1/L\) of the exact finite-volume integral of the sphere [Krüger and Vlugt (2018)].

Parameters:

r (np.ndarray) – Radial distance for weight function.

Returns:

Weight function.

Return type:

np.ndarray

Notes

The weight function, \(u_1(r)\), is defined as:

\[u_1(r) = 4 \pi r^2 \left( 1 - x^3 \right).\]
u2_weight(r: ndarray) → ndarray[source]

Approximation of KBI in the thermodynamic limit using the exact geometric weighting function for hyperspheres and approximation of the surface term [Krüger and Vlugt (2018)].

Parameters:

r (np.ndarray) – Radial distance for weight function.

Returns:

Weight function.

Return type:

np.ndarray

Notes

The weight function, \(u_2(r)\), is defined as:

\[\begin{split}\begin{aligned} u_2(r) &= w(r) \left( 1 + \frac{3}{2}x + \frac{9}{4}x^2 \right), \\ &= 4 \pi r^2 \left( 1 - \frac{23}{8}x^3 + \frac{3}{4}x^4 + \frac{9}{8}x^5 \right). \end{aligned}\end{split}\]
compute_running_kbi(weight_type: str)[source]

Compute KBI as a function of radial distance between molecules, also referred to as the running KBI.

We follow the different correction schemes dicussed in Krüger and Vlugt (2018), where various weight_type implement different corrections.

Parameters:

weight_type (str) – Type of weight function for finite-volume corrections. Options: (‘none’, ‘u0’, ‘u1’, ‘u2’, ‘geometric’)

Returns:

KBI values as a numpy array corresponding to distances \(r\) from the RDF.

Return type:

np.ndarray

Notes

For weight_type='geometric', the KBI is computed using the formula:

\[\begin{split}\begin{aligned} G^V(L) &= \int_0^L h(r) w(r) dr, \\ &= G^\infty + \frac{1}{L}F^\infty + \mathcal{O}(L^{-2}), \end{aligned}\end{split}\]
where:

  • \(L \equiv 6V/A\) is the linear dimension of the system.

  • \(F^\infty\) being the surface term, related to surface effects of the subvolume

  • \(h(r)\) is the total correlation function.

While the above equation provides an exact relation to \(G^\infty\) it requires a linear region for extrapolation to the thermodynamic limit, which can be a disadvantage if a linear region is not identified. Krüger and Vlugt proposed a method for direct estimation of \(G^\infty\), where the accuracy is depending upon weight functions of the form, \(u_k(r)\):

\[G^\infty \approx G_k(L) = \int_0^L h(r) u_k(r) dr.\]
property running_kbi_map: dict[str, ndarray]

Compute the running KBI for each weight_type.

Type:

dict[str, np.ndarray]

property rkbi: ndarray

Running KBI for class initialized weight_type.

Type:

np.ndarray

compute_geometric_extrapolation(positions: tuple[float, float] | None = None, maximize_r2: bool = True, r2_threshold: float = 0.99, min_r_range: float = 1.0, errors: Literal['raise', 'warn', 'ignore'] = 'raise', store: bool = True) → dict[source]

Extrapolate KBI to the thermodynamic limit for weight_type='geometric' using linear regression [Dawass, Krüger, et al. (2020)].

Parameters:

  • position (tuple[float,float], optional) – Range of r values to include for linear fit. Values outside this range will be excluded.

  • maximize_r2 (bool, optional) – Search through valid positions to find the r values that maximize \(R^2\) with a r range greater than 1.0 nm, and that include the last 10% of data (this ensures some range in the beginning is not selected).

  • r2_threshold (float, optional) – Set a \(R^2\) threshold to satisfy KBIConvergence.

  • min_r_range (float, optional) – Minimum range of r values to be required for KBIConvergence (Only for maximize_r2=True).

  • errors (Literal["raise", "warn", "ignore"], optional) – Error mode for handling KBIConvergenceErrors. Only for weight_type='geometric'.

  • store (bool, optional) – Decides if result will be stored in geometric_extrapolation_result.

Returns:

Dictionary containing results from linear extrapolation.

Return type:

dict

Notes

Linear extrapolation of KBI is only performed for the weight_type='geometric' and provides the exact value in the thermodynamic limit. The relation to finite-volume KBI is given through:

\[L G^V(L) = L G^\infty + F^\infty,\]

where convergence is met if a linear region greater than min_r_range is identified and \(R^2 >\) r2_threshold.

property geometric_extrapolation_result: dict

By default, returns the result of the linear extrapolation using maximize_r2=True. This can be updated by running compute_geometric_extrapolation with different arguments and setting store=True.

Type:

dict

compute_kbi(weight_type: str) → float[source]

Get KBI value in the thermodynamic limit according to weight_type.

Parameters:

weight_type (str) – Type of weight function for finite-volume corrections. Options: (‘none’,’u0’,’u1’,’u2’,’geometric’)

Returns:

KBI in the thermodynamic limit

Return type:

float

Notes

For weight_type='geometric'; KBI is compute via extrapolation to the thermodynamic limit. (see also: compute_geometric_extrapolation())

\[L G^V(L) = L G^\infty + F^\infty\]

For all other type of weight_type, KBI is approximated by averaging the value of \(G_k(L)\) over the last 0.1 nm.

\[G^\infty \approx G_k(L)\]
property kbi: float

KBI value in the thermodynamic limit for class initialized weight_type.

Type:

float

add_kbi(axhandle: axes, weight_type: str, **kwargs) → None[source]

Add KBI as a function of radial distance to axes.

Parameters:

  • axhandle (plt.axes) – Axes to add plot.

  • weight_type (str) – Type of weight function for finite-volume corrections. Options: (‘none’,’u0’,’u1’,’u2’,’geometric’)

add_lkbi(axhandle: axes, **kwargs) → None[source]

Add scaled KBI as a function of radial distance to axes, for plotting geometric linear extrapolation.

Parameters:

axhandle (plt.axes) – Axes to add plot.

add_lkbi_fit(axhandle: axes, result: dict | None = None, **kwargs) → None[source]

Add fit results from geometric linear extrapolation to axes.

Parameters:

  • axhandle (plt.axes) – Axes to add plot.

  • result (dict, optional) – Geometric extrapolation result.

add_kbi_value(axhandle: axes, weight_type: str, **kwargs) → None[source]

Add horizontal line for KBI in the thermodynamic limit for a given weight_type.

Parameters:

  • axhandle (plt.axes) – Axes to add plot.

  • weight_type (str) – Type of weight function for finite-volume corrections. Options: (‘none’,’u0’,’u1’,’u2’,’geometric’)

plot_kbi(weight_type: str | None = None, cmap: str | None = None, add_kbi_value: bool = False, filepath: str | None = None, **kwargs) → None[source]

Create figure for plotting single KBI for a given weight type.

Parameters:

  • weight_type (str, optional) – Weight function to plot. Options: (‘none’,’u0’,’u1’,’u2’,’geometric’). Default will go to class initialized weight_type.

  • cmap (str, optional) – Matplotlib colormap name.

  • add_kbi_value (bool, optional) – Add horizontal value for KBI with weight_type='geometric'.

  • filepath (str, optional) – Filepath to save figure as.

plot_kbi_compare(weight_types: list[str], cmap: str | None = None, add_kbi_value: bool = False, filepath: str | None = None) → None[source]

Create figure for comparing KBI’s from various correction methods.

Parameters:

  • weight_types (list[str]) – Weight functions to include in plot. Options: (‘none’,’u0’,’u1’,’u2’,’geometric’)

  • cmap (str, optional) – Matplotlib colormap name.

  • add_kbi_value (bool, optional) – Add horizontal value for KBI with weight_type='geometric'.

  • filepath (str, optional) – Filepath to save figure as.

plot_kbi_extrapolation(result: dict | None = None, color: str = 'turquoise', linewidth: float = 1.0, fitcolor: str = 'k', fitstyle: tuple | str = (0, (3, 2)), fitwidth: float = 2, filepath: str | None = None) → None[source]

Create figure for evaluating geometric extrapolation.

Parameters:

  • result (dict, optional) – Geometric extrapolation result.

  • color (str, optional) – Color for KBI value.

  • linewidth (float, optional) – Linewidth for plotting KBI value.

  • fitcolor (str, optional) – Color for linear fit.

  • fitstyle (tuple | str, optional) – Linestyle for linear fit.

  • fitwidth (float, optional) – Linewidth for linear fit.

  • filepath (str, optional) – Filepath to save figure as.

plot_kbi_compare_extrapolation(weight_types: list[str], cmap: str | None = None, add_kbi_value: bool = False, filepath: str | None = None)[source]

Combine KBI comparisons of weight_type and geometric extrapolation on a single plot.

Parameters:

  • weight_types (list[str]) – Weight functions to include in plot. Options: (‘none’,’u0’,’u1’,’u2’,’geometric’)

  • cmap (str, optional) – Matplotlib colormap name.

  • add_kbi_value (bool, optional) – Add horizontal value for KBI with weight_type='geometric'.

  • filepath (str, optional) – Filepath to save figure as.

References