Kirkwood-Buff Thermodynamics and Structure Factors

Compute thermodynamic properties and structure factors from Kirkwood-Buff integrals (KBIs) across multicomponent systems.

KBThermo applies Kirkwood-Buff theory to a matrix of pairwise KB integrals and constructs thermodynamic property matrices such as:

  • hessians of Gibbs mixing free energy,

  • activity coefficient derivatives,

  • decouples enthalpic vs. entropic contribution to Gibbs mixing free energy,

  • structure factors (partial, Bhatia-Thornton),

  • and related x-ray intensities.

The class operates at constant temperature and uses system metadata (densities, compositions, species identities) provided by a SystemCollection object. It supports multiple strategies for integrating activity coefficient derivatives, including numerical integration and polynomial fitting.

Note

  • KBThermo does not compute KB integrals itself; it consumes a precomputed KBI matrix (e.g., from KBICalculator).

  • All thermodynamic quantities are computed consistently across mixtures, enabling comparison of multicomponent systems or concentration series.

  • Designed for automated workflows within the KBKit analysis pipeline.

class kbkit.kbi.thermodynamics.KBThermo(systems: SystemCollection, kbi: PropertyResult, activity_integration_type: Literal['numerical', 'polynomial'] = 'numerical', activity_polynomial_degree: int = 5)[source]

Bases: object

Apply Kirkwood-Buff (KB) theory to calculate thermodynamic properties.

This class inherits system properties from SystemCollection and uses them for the calculation of thermodynamic properties.

Parameters:

  • systems (SystemCollection) – SystemCollection at a constant temperature.

  • kbi (PropertyResult) – KBI values for each pairwise interaction.

  • activity_integration_type (str, optional) – Method for performing integration of activity coefficient derivatives.

  • activity_polynomial_degree (int, optional) – Polynomial degree for fitting activity coefficient derivatives, if activity_integration_type is polynomial.

kbi(units: str = 'nm^3/molecule') → ndarray[source]

KBI values in desired units.

R(units: str = 'kJ/mol/K') → float[source]

float: Gas constant.

temperature(units: str = 'K') → ndarray[source]

np.ndarray: 1D array of Temperatures of each system.

RT(units: str = 'kJ/mol') → ndarray[source]

np.ndarray: Gas constant (kJ/mol/K) x simulation Temperature.

rho(units: str = 'molecule/nm^3') → ndarray[source]

np.ndarray: 1D array of number density of each system.

v_bar(units: str = 'cm^3/mol') → ndarray[source]

Ideal molar volumes.

\[\bar{V} = \sum_i x_i \bar{V}_i^{pure}\]
Return type:

np.ndarray

property z_i: ndarray

Electrons present in the system mapped to molecules.

Type:

np.ndarray

property z_bar: ndarray

Electrons as a function of composition.

\[\bar{Z} = \sum_i x_i Z_i\]
Return type:

np.ndarray

property z_i_diff: ndarray

Difference in electrons from the last element.

\[\Delta Z_i = Z_i - Z_n\]
where:

  • \(Z_n\): Last element in Z_i()

from \(i=1 \rightarrow n-1\) where \(n\) is the number of molecule types present.

Return type:

np.ndarray

property delta_ij: ndarray

Kronecker delta between pairs of unique molecules (n x n array).

Type:

np.ndarray

B() → ndarray[source]

Calculates the fluctuation matrix, B.

The matrix B represents particle number fluctuations within a fixed volume. In Kirkwood-Buff theory, it serves as the thermodynamic bridge between KBI (G) and the Helmholtz Hessian (A).

Returns:

A 3D matrix with shape (n_sys, n_i, n_i), where n_sys is the number of systems and n_i is the number of unique components.

Return type:

np.ndarray

\[B_{ij} = (A^{-1})_{ij} = \frac{\langle \Delta N_i \Delta N_j \rangle}{V} = x_i \delta_{ij} + \rho x_i x_j G_{ij}\]
where:

  • \(\rho\): Mixture number density.

  • \(x_i\): Mole fraction of species \(i\).

  • \(G_{ij}\): Kirkwood-Buff integral (KBI) for the pair \(i,j\).

  • \(\delta_{ij}\): Kronecker delta (\(\delta_{ij}=1\) if \(i=j\), else 0).

Note

This matrix describes the system in the grand canonical (\(\mu VT\)) limit. It is the direct mathematical inverse of the Helmholtz Hessian at constant volume (\(B = A^{-1}\)).

A() → ndarray[source]

Calculates the Helmholtz Hessian matrix, A.

The Helmholtz Hessian consists of the second partial derivatives of the Helmholtz free energy with respect to particle numbers. It is computed via the inverse of the fluctuation matrix, B.

Returns:

A 3D matrix with shape (n_sys, n_i, n_i), where n_sys is the number of systems and n_i is the number of unique components.

Return type:

np.ndarray

\[A_{ij} = (B^{-1})_{ij} = \frac{1}{RT} \left( \frac{\partial \mu_i}{\partial N_j} \right)_{T,V,N_{k \neq j}}\]
where:

  • \(\mu_i\): Chemical potential of species \(i\).

  • \(N_j\): Particle number (moles) of species \(j\).

Note

This matrix corresponds to the constant volume (canonical) ensemble. Unlike the Gibbs Hessian, this is an $n times n$ matrix and is directly related to particle number fluctuations.

M(units: str = 'kJ/mol') → ndarray[source]

Calculates the full curvature matrix, M.

The matrix M represents the second derivatives of the Gibbs free energy with respect to particle numbers. It is the unconstrained \(n \times n\) Hessian that accounts for volume fluctuations by applying the constant pressure correction to the Helmholtz Hessian.

Returns:

A 3D matrix with shape (n_sys, n_i, n_i), where n_sys is the number of systems and n_i is the number of unique components.

Return type:

np.ndarray

\[\begin{split}\begin{aligned} M_{ij} &= \left(\frac{\partial \mu_i}{\partial N_j}\right)_{T,P,N_{k \neq j}} \\ &= \left(\frac{\partial \mu_i}{\partial N_j}\right)_{T,V,N_{k \neq j}} - \frac{\bar{V}_i \bar{V}_j}{\bar{V} \kappa_T} \\ &= RT \left[ A_{ij} - \frac{\left(\sum_{k=1}^n x_k A_{ik}\right) \left(\sum_{k=1}^n x_k A_{jk}\right)}{\sum_{m=1}^n\sum_{n=1}^n x_m x_n A_{mn}} \right] \\ \end{aligned}\end{split}\]
where:

  • \(A_{ij}\): Elements of the Helmholtz Hessian.

  • \(x_k\): Mole fraction of molecule \(k\).

  • \(\bar{V}_i\): Partial molar volume of molecule \(i\).

  • \(\bar{V}\): Molar volume of the mixture.

  • \(\kappa_T\): Isothermal compressibility of the mixture.

Note

This matrix corresponds to the constant pressure (isobaric-isothermal) ensemble. Due to the Gibbs-Duhem relation, this matrix is mathematically singular (\(\det(M) = 0\)).

isothermal_compressibility(units: str = '1/kPa') → ndarray[source]

Isothermal compressibility, \(\kappa_T\), of mixture.

Returns:

A 1D array with shape (n_sys), where n_sys is the number of systems.

Return type:

np.ndarray

\[\kappa_T RT = \frac{1}{\rho \sum_{j=1}^n \sum_{k=1}^n x_j x_k A_{jk}}\]
where:

  • \(\rho\): Mixture number density.

  • \(A_{ij}\): Element of Helmholtz Hessian matrix for molecules \(i,j\).

molar_volume(units: str = 'cm^3/mol') → ndarray[source]

Partial molar volume of individual components.

Parameters:

units (str, optional) – Desired output units. Defaults to “cm^3/mol”.

Returns:

A 2D array with shape (n_sys, n_i), where n_sys is the number of systems and n_i is the number of unique components.

Return type:

np.ndarray

\[\bar{V}_i = \frac{\sum_{j=1}^n x_j A_{ij}}{\rho \sum_{j=1}^n \sum_{k=1}^n x_j x_k A_{jk}}\]
where:

  • \(\rho\): Mixture number density.

  • \(A_{ij}\): Element of Helmholtz Hessian matrix for molecules \(i,j\).

H(units: str = 'kJ/mol') → ndarray[source]

Calculates the Hessian matrix of Gibbs mixing free energy, H.

Returns:

A 3D matrix with shape (n_sys, n_i-1, n_i-1), where n_sys is the number of systems and n_i is the number of unique components.

Return type:

np.ndarray

\[H_{ij} = \left( \frac{\partial \mu_i}{\partial x_j} \right)_{T,P} = M_{ij} - M_{in} - M_{jn} + M_{nn}\]
where:

  • \(M_{ij}\): Element of the curvature matrix M for molecules \(i,j\)

Note

This matrix is defined in the (n-1) x (n-1) composition space. It represents the stability of the system at constant pressure. A state is considered stable only if this matrix is positive definite.

det_H(units: str = 'kJ/mol') → ndarray[source]

Calculates the determinant of the Gibbs Hessian, H.

The determinant of the reduced Hessian matrix is the primary indicator of thermodynamic stability in a multicomponent mixture.

Returns:

A 1D array of shape (n_sys)

Return type:

np.ndarray

\[\mathcal{D} = \det(\mathbf{H})\]

Note

A system is thermodynamically stable if and only if \(\det(H) > 0\). The condition \(\det(H) = 0\) defines the spinodal line, beyond which the single-phase mixture becomes spontaneously unstable to infinitesimal fluctuations.

dmu_dxi(units: str = 'kJ/mol') → ndarray[source]

Calculates the chemical potential derivatives with respect to mole fraction.

This property returns the diagonal elements of the chemical potential derivative matrix, representing the response of each species’ chemical potential to its own mole fraction, constrained by the Gibbs-Duhem relation.

Returns:

A 2D array of shape (n_sys, n_i), where n_sys is the number of systems and n_i is the number of unique components.

Return type:

np.ndarray

\[\left(\frac{\partial \mu_i}{\partial x_i}\right)_{T,P} = \sum_j \left(\frac{\partial \mu_i}{\partial n_i}\right)_{T,P} \delta_{ij}\]

Note

For the first \(n-1\) components, the derivatives are transformed from the particle-number basis (M). The \(n^{th}\) component is calculated via the Gibbs-Duhem equation:

\[\frac{\partial \mu_n}{\partial x_n} = \frac{1}{x_n} \sum_{j=1}^{n-1} x_j \frac{\partial \mu_j}{\partial x_j}\]

This ensures that the composition derivatives are thermodynamically consistent across the entire mixture.

Gamma() → ndarray[source]

Calculates the thermodynamic factors, \(\Gamma_i\), demonstrated by Fingerhut et al. (2020).

The thermodynamic factor scales the composition dependence of the chemical potential.

Returns:

A 3D array of shape (n_sys, n_i, n_i).

Return type:

np.ndarray

\[\Gamma_{ij} = \frac{x_i}{RT} \left( \frac{\partial \mu_i}{\partial x_j} \right)_{T,P}\]
dlngamma_dxi() → ndarray[source]

Derivative of natural logarithm of the activity coefficient of molecule \(i\) with respect to its own mole fraction.

This represents the deviation from ideality in the chemical potential gradient.

Returns:

A 3D matrix with shape (n_sys, n_i, n_i)

Return type:

np.ndarray

\[\frac{\partial \ln{\gamma_i}}{\partial x_i} = \frac{1}{R T}\left(\frac{\partial \mu_i}{\partial x_i}\right)_{T,P} - \frac{1}{x_i}\]
where:

  • \(\mu_i\): Chemical potential of molecule \(i\)

  • \(\gamma_i\): Activity coefficient of molecule \(i\)

  • \(x_i\): Mole fraction of molecule \(i\)

lngamma() → ndarray[source]

Natural logarithm of activity coefficients.

Integrate the derivative of activity coefficients to obtain \(\ln{\gamma_i}\) for each component. Use either numerical methods (trapezoidal rule) or polynomial fitting for integration. These parameters are chosen by the activity_integration_type and activity_polynomial_degree in KBThermo initialization.

Returns:

A 2D array with shape (n_sys, n_i)

Return type:

np.ndarray

Notes

The general formula for activity coefficient integration is:

\[\ln{\gamma_i}(x_i) = \int \left(\frac{\partial \ln{\gamma_i}}{\partial x_i}\right) dx_i\]

Polynomial integration: the method fits a polynomial, \(P(x_i)\), to the derivative data and integrates:

\[\ln{\gamma_i}(x_i) = \int P(x_i) dx_i\]

The integration constant is chosen so that \(\ln{\gamma_i}\) obeys the boundary condition at the reference state.

Numerical Integration: The trapezoidal rule is used to approximate the integral because an analytical solution is not available. The integral is approximated as:

\[\ln{\gamma_i}(x_i) \approx \sum_{a=a_0}^{N-1} \frac{(x_i)_{a+1}-(x_i)_a}{2} \left[\left(\frac{\partial \ln{\gamma_i}}{\partial x_i}\right)_{a} + \left(\frac{\partial \ln{\gamma_i}}{\partial x_i}\right)_{a+1}\right]\]
where:

  • \(\ln{\gamma_i}(x_i)\): Natural logarithm of the activity coefficient of molecule i at mole fraction \(x_i\).

  • \(a\): Index of summation.

  • \(a_0\): Starting value for index of summation.

  • \(N\): Number of data points to sum over.

  • \(x_i\): Mole fraction of component \(i\).

  • \(\left(\frac{\partial \ln{\gamma_i}}{\partial x_i}\right)_{a}\): Derivative of the natural logarithm of the activity coefficient of component i with respect to its mole fraction, evaluated at point a.

The integration starts at a reference state where \(x_i = a_0\) and \(\ln{\gamma_i}(a_0) = 0\).

property activity_metadata: ActivityMetadata

Container for results from activity coefficient integration.

Type:

ActivityMetadata

g_ex(units: str = 'kJ/mol') → ndarray[source]

Gibbs excess energy from activity coefficients.

\[\frac{G^E}{RT} = \sum_{i=1}^n x_i \ln{\gamma_i}\]
where:

  • \(x_i\): Mole fraction of molecule \(i\).

  • \(\gamma_i\): Activity coefficient of molecule \(i\).

h_mix(units: str = 'kJ/mol') → ndarray[source]

Enthalpy of mixing. Requires pure component simulations.

\[\Delta H_{mix} = H - \sum_{i} x_i H_i^{pure}\]
where:

  • \(H\): Enthalpy directly from simulation.

  • \(H_i^{pure}\): Enthalpy directly from simulation for pure \(i\).

s_ex(units: str = 'kJ/mol/K') → ndarray[source]

Excess entropy from mixing enthalpy and Gibbs excess energy. Requires pure component simulations.

\[S^E = \frac{\Delta H_{mix} - G^E}{T}\]
where:

  • \(x_i\): Mole fraction of molecule \(i\).

  • \(\Delta H_{mix}\): Enthalpy of mixing.

  • \(G^E\): Excess Gibbs energy.

g_id(units: str = 'kJ/mol') → ndarray[source]

Ideal free energy calculated from mole fractions.

\[\frac{G^{id}}{RT} = \sum_{i=1}^n x_i \ln{x_i}\]
where:

  • \(x_i\) is mole fraction of molecule \(i\).

s_mix(units: str = 'kJ/mol/K') → ndarray[source]

Mixing entropy, requires pure component simulations.

\[\begin{split}\begin{aligned} \Delta S_{mix} &= S^E + S^{id} \\ &= S^E - R \sum_{i=1}^n x_i \ln{x_i} \end{aligned}\end{split}\]
where:

  • \(x_i\): Mole fraction of molecule \(i\).

  • \(S^E\): Excess entropy.

  • \(S^{id}\): Ideal entropy.

g_mix(units: str = 'kJ/mol') → ndarray[source]

Gibbs mixing free energy calculated from excess and ideal contributions.

\[\begin{split}\begin{aligned} \Delta G_{mix} &= G^E + G^{id} \\ &= \Delta H_{mix} - T \Delta S_{mix} \end{aligned}\end{split}\]
where:

  • \(x_i\): Mole fraction of molecule \(i\).

  • \(G^E\): Excess Gibbs energy.

  • \(\Delta G_{mix}\): Gibbs free energy of mixing.

  • \(\Delta H_{mix}\): Enthalpy of mixing.

  • \(\Delta S_{mix}\): Entropy of mixing.

s0_ij() → ndarray[source]

Partial structure factors for pairwise interaction between components.

Notes

Partial structure factor, \(\hat{S}_{ij}(0)\), is calculated via:

\[\hat{S}_{ij}(0) = B_{ij} = \rho x_i x_j G_{ij} + x_i \delta_{i,j}\]
where:

  • \(G_{ij}\) is the KBI value for molecules \(i,j\).

  • \(B\) is the fluctuation matrix.

Note

Note that the normalization used here differs from that of the Ashcroft-Langreth partial structure factors used in some texts.

s0_x() → ndarray[source]

Bhatia-Thornton composition-composition structure factor as q \(\rightarrow\) 0, extended to a multicomponent system.

Notes

Structure factor, \(\hat{S}_{ij}^{x}(0)\), is a 3D matrix (composition x n-1 x n-1 components) and is calculated via:

\[\hat{S}_{ij}^{x}(0) = \hat{S}_{ij}(0) - x_i \sum_{k=1}^n \hat{S}_{kj}(0) - x_j \sum_{k=1}^n \hat{S}_{ki}(0) + x_i x_j \sum_{k=1}^n \sum_{l=1}^n \hat{S}_{kl}(0)\]

for i and j from 1 to n-1.

s0_xp() → ndarray[source]

Bhatia-Thornton composition-density structure factor as q \(\rightarrow\) 0, extended to a multicomponent system.

Notes

Structure factor, \(\hat{S}_{i}^{x\rho}(0)\), is a 2D array (composition x n-1 components) and is calculated via:

\[\hat{S}_{i}^{x\rho}(0) = \sum_{k=1}^n \hat{S}_{ik}(0) - x_i \sum_{k=1}^n \sum_{l=1}^n \hat{S}_{kl}(0)\]

for i from 1 to n-1.

s0_p() → ndarray[source]

Bhatia-Thornton density-density structure factor as q \(\rightarrow\) 0, extended to a multicomponent system.

Notes

Structure factor, \(\hat{S}^{\rho}(0)\), is a 1D vector (composition) and is calculated via:

\[\hat{S}^{\rho}(0) = \sum_{k=1}^n \sum_{l=1}^n \hat{S}_{kl}(0)\]
s0_kappa() → ndarray[source]

Contribution from isothermal compressibility to Bhatia-Thornton density-density structure factor as q \(\rightarrow\) 0.

Notes

Structure factor, \(\hat{S}^{\kappa_T}(0)\), is calculated via:

\[\hat{S}^{\kappa_T}(0) = \frac{RT \kappa_T}{\bar{V}}\]
s0_x_e() → ndarray[source]

Contribution from extended Bhatia-Thornton composition-composition structure factors to electron density structure factor as q \(\rightarrow\) 0.

Notes

Structure factor, \(\hat{S}^{x,e}(0)\), is a 1D vector (composition) and is calculated via:

\[\hat{S}^{x,e}(0) = \sum_{i=1}^{n-1}\sum_{j=1}^{n-1} \left( Z_i - Z_n \right) \left( Z_j - Z_n \right) \hat{S}_{ij}^{x}(0)\]
s0_xp_e() → ndarray[source]

Contribution from extended Bhatia-Thornton composition-density structure factors to electron density structure factor as q \(\rightarrow\) 0.

Notes

Structure factor, \(\hat{S}^{x\rho,e}(0)\), is a 1D vector (composition) and is calculated via:

\[\hat{S}^{x\rho,e}(0) = 2 \bar{Z} \sum_{i=1}^{n-1} \left( Z_i - Z_n \right) \hat{S}_{i}^{x\rho}(0)\]
s0_p_e() → ndarray[source]

Contribution from extended Bhatia-Thornton density-density structure factors to electron density structure factor as q \(\rightarrow\) 0.

Notes

Structure factor, \(\hat{S}^{\rho,e}(0)\), is a 1D vector (composition) and is calculated via:

\[\hat{S}^{\rho,e}(0) = \bar{Z}^2 \hat{S}^{\rho}(0)\]
s0_kappa_e() → ndarray[source]

Contribution from isothermal compressibility part of Bhatia-Thornton density-density structure factor to electron density structure factor as q \(\rightarrow\) 0.

Notes

Structure factor, \(\hat{S}^{\kappa_T, e}(0)\), is calculated via:

\[\hat{S}^{\kappa_T, e}(0) = \bar{Z}^2 \hat{S}^{\kappa_T}(0)\]
s0_e() → ndarray[source]

Electron density structure factor as q \(\rightarrow\) 0 for the entire mixture.

Notes

Structure factor, \(\hat{S}^{e}(0)\), can be calculated via partial or from Bhatia-Thornton structure factors (both are equivalent):

\[\begin{split}\begin{aligned} \hat{S}^{e}(0) &= \sum_{i=1}^n \sum_{j=1}^n Z_i Z_j \hat{S}_{ij}(0) \\ &= \hat{S}^{x,e}(0) + \hat{S}^{x\rho,e}(0) + \hat{S}^{\rho,e}(0) \end{aligned}\end{split}\]
i0_x(units: str = '1/cm') → ndarray[source]

Contribution from extended Bhatia-Thornton composition-composition structure factors to x-ray intensity as q \(\rightarrow\) 0.

Notes

X-ray intensity, \(I^{x}(0)\), is calculated via:

\[I^{x}(0) = r_e^2 \rho N_A \hat{S}^{x,e}(0)\]
i0_xp(units: str = '1/cm') → ndarray[source]

Contribution from extended Bhatia-Thornton composition-density structure factors to x-ray intensity as q \(\rightarrow\) 0.

Notes

X-ray intensity, \(I^{x\rho}(0)\), is calculated via:

\[I^{x\rho}(0) = r_e^2 \rho N_A \hat{S}^{x\rho,e}(0)\]
i0_p(units: str = '1/cm') → ndarray[source]

Contribution from extended Bhatia-Thornton density-density structure factors to x-ray intensity as q \(\rightarrow\) 0.

Notes

X-ray intensity, \(I^{\rho}(0)\), is calculated via:

\[I^{\rho}(0) = r_e^2 \rho N_A \hat{S}^{\rho,e}(0)\]
i0_kappa(units: str = '1/cm') → ndarray[source]

Contribution from isothermal compressibility part of Bhatia-Thornton density-density structure factor to x-ray intensity as q \(\rightarrow\) 0.

Notes

X-ray intensity, \(I^{\kappa_T}(0)\), is calculated via:

\[I^{\kappa_T}(0) = r_e^2 \rho N_A \hat{S}^{\kappa_T,e}(0)\]
i0(units: str = '1/cm') → ndarray[source]

X-ray intensity as q \(\rightarrow\) 0 for entire mixture.

Notes

X-ray intensity, \(I(0)\), is calculated via:

\[\begin{split}\begin{aligned} I(0) &= r_e^2 \rho N_A \hat{S}^e \\ &= I^x(0) + I^{x\rho}(0) + I^{\rho}(0) \end{aligned}\end{split}\]
property results: dict[str, PropertyResult]

Container for PropertyResult objects for KBI and KBI-derived quantities.

Type:

dict

plotter(molecule_map: dict[str, str] | None = None) → ThermoPlotter[source]

Create a ThermoPlotter for visualizing KBI and KBI-derived properties as a function of composition.

Returns:

Plotter instance for computing KBI-derived thermodynamic properties.

Return type:

ThermoPlotter

References