chemistrykit.statmech#

chemistrykit.statmech hub · History · Examples · API reference

chemistrykit.statmech: statistical mechanics of molecules.

Translational/rotational/vibrational partition functions and the thermodynamic functions (U, S, Cv) derived from them (including the Sackur-Tetrode translational entropy); the Maxwell-Boltzmann speed distribution; a canonical-ensemble lattice-gas adsorption model (a statistical-mechanical derivation of the Langmuir isotherm); the Debye heat capacity of a solid; exact Ising-model results (the 1D transfer matrix, Kramers-Wannier duality, Onsager’s 2D solution); and the second virial coefficient from a pair potential.

chemistrykit.md’s trajectories are a natural cross-check for this domain’s distribution functions – see examples/md/lj_fluid/plot_02_maxwell_boltzmann_check.py.

class chemistrykit.statmech.DebyeSolid(debye_temperature)[source]#

Bases: object

Debye’s continuum model of a monatomic crystal’s lattice vibrations.

Parameters:: debye_temperature (float) – The Debye temperature \(\Theta_D=h\nu_D/k_B\), in K.

Examples

Far above \(\Theta_D\) the molar heat capacity approaches the Dulong-Petit value \(3R\); far below, it follows Debye’s \(T^3\) law:

>>> copper = DebyeSolid(debye_temperature=343.0)
>>> round(float(copper.heat_capacity_v(5000.0) / (3 * 8.31446261815324)), 3)
1.0
>>> ratio = copper.heat_capacity_v(3.0) / copper.low_temperature_heat_capacity(3.0)
>>> round(float(ratio), 6)
1.0

heat_capacity_v(T, N=6.02214076e+23)[source]#

Return the constant-volume heat capacity of N atoms (default one mole).

Parameters:

T (float or array_like) – Absolute temperature, in K (must be positive).
N (float) – Number of atoms.

Returns:

float or numpy.ndarray – Heat capacity, in J/K.

internal_energy(T, N=6.02214076e+23)[source]#

Return the thermal internal energy \(U(T)-U(0)\) of N atoms.

\(U-U_0 = 9Nk_BT(T/\Theta_D)^3\int_0^{\Theta_D/T}x^3/(e^x-1)\,dx\), measured from the zero-point energy \(U_0=\frac98Nk_B\Theta_D\).

Parameters:

T (float or array_like) – Absolute temperature, in K (must be positive).
N (float)

Returns:

float or numpy.ndarray – Energy, in J.

low_temperature_heat_capacity(T, N=6.02214076e+23)[source]#

Return Debye’s \(T^3\) law, \(C_V=\frac{12\pi^4}{5}Nk_B(T/\Theta_D)^3\), exact as \(T\to0\).

Parameters:

T (float or array_like) – Absolute temperature, in K.
N (float)

Returns:

float or numpy.ndarray – Heat capacity, in J/K.

class chemistrykit.statmech.IdealGasMolecule(mass, volume, moment_of_inertia=None, symmetry_number=1, vibrational_frequencies=())[source]#

Bases: object

Combined translational + rotational + vibrational thermodynamic functions for one ideal-gas species.

Sums the independent contributions of a TranslationalPartitionFunction, an optional RotationalPartitionFunctionLinear, and zero or more VibrationalPartitionFunctionHarmonic modes – valid because the total molecular partition function factorizes, \(q=q_{\text{trans}}q_{\text{rot}}q_{\text{vib},1}q_{\text{vib},2}\cdots\), whenever these degrees of freedom are independent (the standard approximation for a rigid-rotor/harmonic-oscillator molecule; Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 13.1).

Parameters:

mass (float) – Molecular mass, in kg.
volume (float) – Container volume, in m^3.
moment_of_inertia (float) – Moment of inertia for a linear-rotor rotational contribution; omit for an atom (no rotational degrees of freedom).
symmetry_number (int)
vibrational_frequencies (Sequence[float]) – Vibrational mode frequencies, in Hz.

Examples

At very high temperature every mode approaches its classical equipartition heat capacity: \(\frac32R\) (translation) + \(R\) (linear rotation) + \(R\) per vibrational mode:

>>> molecule = IdealGasMolecule(mass=6.63e-26, volume=1.0e-3, moment_of_inertia=1.45e-46, vibrational_frequencies=[8.7e13])
>>> Cv_over_R = molecule.heat_capacity_v(1.0e7) / 8.31446261815324
>>> round(float(Cv_over_R), 2)
3.5

entropy(T, N=6.02214076e+23)[source]#

Total entropy: sum of each mode’s contribution.

Parameters:: N (float)

heat_capacity_v(T, N=6.02214076e+23)[source]#

Total heat capacity: sum of each mode’s contribution.

Parameters:: N (float)

internal_energy(T, N=6.02214076e+23)[source]#

Total internal energy: sum of each mode’s contribution. See chemistrykit.statmech.core.base_system.PartitionFunction.internal_energy().

Parameters:: N (float)

thermodynamic_functions(T, N=6.02214076e+23)[source]#

Bundle q, U, S, Cv, and A at temperature T into a ThermodynamicFunctions.

Parameters:

T (float)
N (float)

Return type:

ThermodynamicFunctions

Returns:

chemistrykit.statmech.core.base_system.ThermodynamicFunctions

class chemistrykit.statmech.Ising1D(coupling, field=0.0)[source]#

Bases: object

The one-dimensional Ising chain, solved exactly by the transfer matrix.

With \(K=J/k_BT\) and \(b=h/k_BT\), the transfer matrix \(T_{ss'}=e^{Kss'+b(s+s')/2}\) has eigenvalues \(\lambda_\pm=e^K[\cosh b\pm\sqrt{\sinh^2b+e^{-4K}}]\), so a periodic chain of N spins has \(Z_N=\lambda_+^N+\lambda_-^N\).

Parameters:

coupling (float) – Exchange coupling J, in J (positive = ferromagnetic).
field (float) – External field h, in J.

Examples

In zero field the magnetization vanishes at every temperature – the one-dimensional chain never orders:

>>> chain = Ising1D(coupling=1.380649e-23 * 100.0)
>>> float(chain.magnetization(50.0))
0.0

but the correlation length grows as \(\xi\approx e^{2K}/2\) on cooling:

>>> round(float(chain.correlation_length(50.0)), 2)
27.3

correlation_length(T)[source]#

Return the zero-field correlation length in lattice spacings, \(\xi=-1/\ln\tanh K\).

Uses the zero-field ratio \(\lambda_-/\lambda_+=\tanh K\) (the stored field is ignored).

free_energy_per_spin(T)[source]#: Return the thermodynamic-limit free energy per spin, \(f=-k_BT\ln\lambda_+\), in J.

heat_capacity_per_spin(T)[source]#: Return the zero-field heat capacity per spin, \(c=k_BK^2\operatorname{sech}^2K\), in J/K.

magnetization(T)[source]#: Return the magnetization per spin, \(m=\sinh b/\sqrt{\sinh^2b+e^{-4K}}\).

partition_function(T, N)[source]#

Return \(Z_N=\lambda_+^N+\lambda_-^N\) for a periodic ring of N spins.

Parameters:

T (float) – Absolute temperature, in K.
N (int) – Number of spins.

Returns:

float

transfer_matrix_eigenvalues(T)[source]#

Return \((\lambda_+,\lambda_-)\) of the 2x2 transfer matrix at temperature T.

Parameters:: T (float) – Absolute temperature, in K.
Returns:: tuple of float

class chemistrykit.statmech.Ising2DOnsager(coupling)[source]#

Bases: object

Onsager’s exact solution of the zero-field square-lattice Ising model.

With \(K=J/k_BT\) and \(\kappa=2\sinh2K/\cosh^22K\):

\[-\beta f = \ln(2\cosh2K)+\frac1\pi\int_0^{\pi/2} \ln\tfrac12\left(1+\sqrt{1-\kappa^2\sin^2\phi}\right)d\phi\]

and, for \(T<T_c\), Yang’s spontaneous magnetization \(m=[1-\sinh^{-4}2K]^{1/8}\).

Parameters:: coupling (float) – Exchange coupling J, in J (positive).

Examples

>>> model = Ising2DOnsager(coupling=1.380649e-23)  # J/k_B = 1 K
>>> round(model.critical_temperature, 4)
2.2692
>>> round(float(model.spontaneous_magnetization(1.0)), 4)
0.9993
>>> float(model.spontaneous_magnetization(3.0))
0.0
>>> round(float(model.internal_energy_per_spin(model.critical_temperature) / model.coupling), 6)
-1.414214

property critical_temperature: float#: The exact critical temperature, in K. See ising_2d_critical_temperature().

free_energy_per_spin(T)[source]#: Return the free energy per spin, in J (Onsager’s single-integral form).

heat_capacity_per_spin(T)[source]#

Return the heat capacity per spin, in J/K (diverges logarithmically at \(T_c\)).

\[\frac{c}{k_B}=\frac2\pi(K\coth2K)^2\Big\{2K_1(\kappa)-2E_1(\kappa) -(1-\kappa')\big[\tfrac\pi2+\kappa'K_1(\kappa)\big]\Big\}\]

with \(\kappa'=2\tanh^22K-1\) and \(E_1\) the complete elliptic integral of the second kind.

internal_energy_per_spin(T)[source]#

Return \(u=-J\coth2K\,[1+\tfrac2\pi(2\tanh^22K-1)K_1(\kappa)]\), in J.

\(K_1\) is the complete elliptic integral of the first kind; its logarithmic divergence at \(T_c\) is multiplied by zero, so u is continuous there (\(u(T_c)=-\sqrt2J\)).

spontaneous_magnetization(T)[source]#: Return Yang’s spontaneous magnetization per spin, \([1-\sinh^{-4}2K]^{1/8}\) below \(T_c\), else 0.

class chemistrykit.statmech.LatticeGasAdsorption(adsorption_energy, mass, T)[source]#

Bases: object

Non-interacting lattice-gas model of adsorption on M independent sites.

Each site’s single-site grand partition function is \(\xi=1+\lambda e^{\beta\epsilon}\) (empty, or occupied with Boltzmann factor \(e^{\beta\epsilon}\) weighted by the gas-phase fugacity \(\lambda=e^{\beta\mu}\)); since the M sites are independent, \(\Xi=\xi^M\) and the average coverage is \(\theta=\langle N\rangle/M=1/(1+e^{-\beta(\mu+\epsilon)})\). Using the ideal-gas chemical potential \(\mu=k_BT\ln(P\Lambda^3/k_BT)\) (with \(\Lambda\) the thermal de Broglie wavelength, chemistrykit.statmech.utils.thermal_wavelength.thermal_de_broglie_wavelength()) turns this into exactly the Langmuir isotherm \(\theta=P/(P+P_0)\) (see coverage(), p_half()).

Parameters:

adsorption_energy (float) – Binding energy per adsorbed molecule, in J (positive = favorable adsorption).
mass (float) – Mass of the adsorbing gas molecule, in kg.
T (float) – Absolute temperature, in K.

static canonical_entropy(N, M)[source]#

Canonical (fixed-N) combinatorial entropy \(S=k_B\ln\binom{M}{N}\).

The exact configurational entropy of placing N indistinguishable adsorbed molecules on M distinguishable sites with no interactions, computed via the numerically exact chemistrykit.statmech.utils.combinatorics.ln_binomial() (rather than Stirling’s approximation).

Parameters:

N (float) – Number of occupied sites (0 <= N <= M).
M (float) – Total number of sites.

Return type:

float

Returns:

float – Entropy, in J/K.

Examples

Zero at the two fully-ordered extremes (all sites empty, or all full – only one microstate each), and maximal exactly at half filling, where it approaches the large-M Stirling estimate \(Mk_B\ln2\) (McQuarrie, Statistical Mechanics, Ch. 8.1):

>>> LatticeGasAdsorption.canonical_entropy(0, 100)
0.0
>>> LatticeGasAdsorption.canonical_entropy(100, 100)
0.0
>>> import numpy as np
>>> M = 100_000
>>> S_half = LatticeGasAdsorption.canonical_entropy(M // 2, M)
>>> relative_error = abs(S_half - M * 1.380649e-23 * np.log(2)) / (M * 1.380649e-23 * np.log(2))
>>> bool(relative_error < 1e-4)
True

coverage(P)[source]#

Fractional coverage \(\theta(P) = P/(P+P_0)\) – the Langmuir isotherm.

Parameters:: P (float or array-like of float) – Gas pressure, in Pa.
Returns:: float or ndarray

Examples

Coverage is exactly one-half at \(P=P_0\) by construction, and approaches the physical limits of 0 (vacuum) and 1 (saturation) at the extremes:

>>> model = LatticeGasAdsorption(adsorption_energy=3.0e-20, mass=4.65e-26, T=300.0)
>>> round(float(model.coverage(model.p_half())), 10)
0.5
>>> float(model.coverage(0.0))
0.0
>>> bool(model.coverage(1.0e15) > 0.9999)
True

p_half()[source]#

Pressure at half coverage, \(P_0=(k_BT/\Lambda^3)e^{-\epsilon/k_BT}\).

Return type:: float
Returns:: float – Pressure, in Pa.

thermal_wavelength()[source]#

The gas molecule’s thermal de Broglie wavelength at T. See chemistrykit.statmech.utils.thermal_wavelength.thermal_de_broglie_wavelength().

Return type:: float

class chemistrykit.statmech.MaxwellBoltzmannSpeedDistribution(mass, temperature, k_b=1.380649e-23)[source]#

Bases: object

The Maxwell-Boltzmann speed distribution for a gas of mass mass at temperature.

Parameters:

mass (float) – Molecular mass (in kg for SI/real units; in reduced units, e.g. those of chemistrykit.md.systems.lj_fluid.LJFluid, pass the reduced mass and k_b=1.0).
temperature (float) – Absolute temperature.
k_b (float) – Boltzmann constant, in the same unit system as mass/temperature.

Examples

The three characteristic speeds stand in the fixed ratio \(v_p:\bar v:v_{\text{rms}} = \sqrt2:\sqrt{8/\pi}:\sqrt3\), independent of mass or temperature – a purely geometric property of the distribution’s shape (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 1.2b):

>>> import numpy as np
>>> dist = MaxwellBoltzmannSpeedDistribution(mass=6.63e-26, temperature=298.15)
>>> vp, vbar, vrms = dist.most_probable_speed(), dist.mean_speed(), dist.rms_speed()
>>> np.allclose([vbar / vp, vrms / vp], [np.sqrt(8.0 / np.pi) / np.sqrt(2.0), np.sqrt(3.0) / np.sqrt(2.0)])
True

The pdf integrates to 1 over all speeds:

>>> from scipy.integrate import quad
>>> total, _ = quad(dist.pdf, 0.0, np.inf)
>>> round(total, 6)
1.0

cdf(v)[source]#

Cumulative probability \(P(V\leq v)\).

\[F(v) = \mathrm{erf}(av) - \frac{2}{\sqrt\pi}av\,e^{-(av)^2}, \qquad a=\sqrt{\frac{m}{2k_BT}}\]

Parameters:: v (float or array-like of float)
Returns:: float or ndarray

Examples

>>> dist = MaxwellBoltzmannSpeedDistribution(mass=6.63e-26, temperature=298.15)
>>> float(dist.cdf(0.0))
0.0

mean_speed()[source]#

See mean_speed().

Return type:: float

most_probable_speed()[source]#

See most_probable_speed().

Return type:: float

pdf(v)[source]#

Probability density f(v) at speed(s) v.

Parameters:: v (float or array-like of float) – Must be non-negative.
Returns:: float or ndarray

rms_speed()[source]#

See rms_speed().

Return type:: float

sample(n, rng=None)[source]#

Draw n random speeds from this distribution.

Each of the 3 Cartesian velocity components of a classical ideal gas is independently Gaussian-distributed with variance \(k_BT/m\); the speed (the norm of the velocity vector) of such a vector is exactly Maxwell-Boltzmann distributed, which is the simplest way to sample it (no inverse-CDF / rejection-sampling machinery needed).

Parameters:

n (int)
rng (int, numpy.random.Generator, or None)

Returns:

ndarray, shape (n,)

Examples

>>> dist = MaxwellBoltzmannSpeedDistribution(mass=6.63e-26, temperature=298.15)
>>> speeds = dist.sample(200_000, rng=0)
>>> bool(abs(speeds.mean() - dist.mean_speed()) / dist.mean_speed() < 0.01)
True

class chemistrykit.statmech.PartitionFunction[source]#

Bases: ABC

Common interface for a single-mode molecular partition function.

Every method accepts N, the number of independent, non-interacting copies of this mode contributing to an extensive thermodynamic quantity – N moles’ worth of molecules if N is left at its default of Avogadro’s number (giving molar quantities directly), a single molecule’s contribution if N=1, or any other ensemble size.

abstractmethod entropy(T, N=6.02214076e+23)[source]#

Return the entropy S(T) contributed by this mode.

Parameters:

T (float)
N (float)

Returns:

float – Entropy, in J/K.

abstractmethod heat_capacity_v(T, N=6.02214076e+23)[source]#

Return the constant-volume heat capacity Cv(T) contributed by this mode.

Parameters:

T (float)
N (float)

Returns:

float – Heat capacity, in J/K.

helmholtz_free_energy(T, N=6.02214076e+23)[source]#

Return \(A=U-TS\) for this mode, from internal_energy() and entropy().

Parameters:

T (float)
N (float)

Returns:

float – Helmholtz free energy, in J.

abstractmethod internal_energy(T, N=6.02214076e+23)[source]#

Return \(U(T)-U(0) = Nk_BT^2(\partial\ln q/\partial T)_V\).

Parameters:

T (float) – Absolute temperature, in K.
N (float) – Number of molecules/modes.

Returns:

float – Energy, in J (relative to this mode’s own ground state).

abstractmethod value(T)[source]#

Return the partition function q(T) for a single molecule/mode.

Parameters:: T (float) – Absolute temperature, in K.
Returns:: float

class chemistrykit.statmech.RotationalPartitionFunctionLinear(moment_of_inertia, symmetry_number=1)[source]#

Bases: PartitionFunction

The rotational partition function of a linear rigid rotor, high-temperature (classical) limit.

\[q_{\text{rot}} = \frac{T}{\sigma\Theta_{\text{rot}}}, \qquad \Theta_{\text{rot}} = \frac{h^2}{8\pi^2Ik_B}\]

where I is the moment of inertia and \(\sigma\) the rotational symmetry number (2 for a homonuclear diatomic, 1 for heteronuclear). This is the classical (high-temperature) approximation to the exact sum over quantized rotational levels, valid whenever \(T\gg\Theta_{\text{rot}}\) – true at room temperature for essentially every molecule except H2 and its isotopologues (McQuarrie, Statistical Mechanics, Ch. 6.3; Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 13.2c). Because this classical limit is used, internal_energy() and heat_capacity_v() are the exact equipartition values (\(Nk_BT\) and \(Nk_B\) respectively, from 2 rotational quadratic degrees of freedom) rather than a sum over levels – this is the approximation being flagged.

Parameters:

moment_of_inertia (float) – Moment of inertia, in kg m^2.
symmetry_number (int) – Rotational symmetry number \(\sigma\).

Examples

>>> q = RotationalPartitionFunctionLinear(moment_of_inertia=1.45e-46, symmetry_number=1)  # ~HCl-like
>>> round(q.heat_capacity_v(298.15, N=1.0) / 1.380649e-23, 6)
1.0

entropy(T, N=6.02214076e+23)[source]#

Return the entropy S(T) contributed by this mode.

Parameters:

T (float)
N (float)

Returns:

float – Entropy, in J/K.

heat_capacity_v(T, N=6.02214076e+23)[source]#

Return the constant-volume heat capacity Cv(T) contributed by this mode.

Parameters:

T (float)
N (float)

Returns:

float – Heat capacity, in J/K.

internal_energy(T, N=6.02214076e+23)[source]#

Return \(U(T)-U(0) = Nk_BT^2(\partial\ln q/\partial T)_V\).

Parameters:

T (float) – Absolute temperature, in K.
N (float) – Number of molecules/modes.

Returns:

float – Energy, in J (relative to this mode’s own ground state).

property rotational_temperature: float#: The rotational temperature \(\Theta_{\text{rot}}=h^2/(8\pi^2Ik_B)\), in K.

value(T)[source]#

Return the partition function q(T) for a single molecule/mode.

Parameters:: T (float) – Absolute temperature, in K.
Returns:: float

class chemistrykit.statmech.ThermodynamicFunctions(T, q, U, S, Cv, A)[source]#

Bases: object

A snapshot of one mode’s (or a combined molecule’s) thermodynamic functions at a given T.

Returned by chemistrykit.statmech.systems.partition_functions.IdealGasMolecule.thermodynamic_functions().

Parameters:

T (float)
q (float)
U (float)
S (float)
Cv (float)
A (float)

A: float#

Helmholtz free energy, U - T*S, in J.

Type:: float

Cv: float#

Constant-volume heat capacity, in J/K.

Type:: float

S: float#

Entropy, in J/K.

Type:: float

T: float#

Temperature, in K.

Type:: float

U: float#

Internal energy (relative to the ground state), in J.

Type:: float

q: float#

The (combined) partition function value at T.

Type:: float

class chemistrykit.statmech.TranslationalPartitionFunction(mass, volume)[source]#

Bases: PartitionFunction

The translational partition function of a particle in a box (volume) V.

\[q_{\text{trans}} = \left(\frac{2\pi mk_BT}{h^2}\right)^{3/2}V\]

the high-temperature (particle-in-a-box energy levels essentially continuous) limit of the exact 3D particle-in-a-box partition function, valid whenever \(V^{1/3}\gg\Lambda\) (the thermal de Broglie wavelength, chemistrykit.statmech.utils.thermal_wavelength.thermal_de_broglie_wavelength()) – true for any real gas well above its condensation point (McQuarrie, Statistical Mechanics, Ch. 6.1). Because translational motion permutes identical molecules through the same physical states (unlike the rotational/vibrational modes below, which are internal to each molecule), the N-molecule thermodynamic functions here include the \(1/N!\) indistinguishability correction, giving the Sackur-Tetrode entropy (see entropy(), sackur_tetrode_entropy()).

Parameters:

mass (float) – Molecular mass, in kg.
volume (float) – Container volume, in m^3.

Examples

Internal energy is exactly the equipartition value \(\frac32Nk_BT\) (3 translational quadratic degrees of freedom) at every temperature – the classical limit is essentially exact here:

>>> q = TranslationalPartitionFunction(mass=6.63e-26, volume=1.0e-3)
>>> round(q.internal_energy(298.15, N=1.0) / (1.5 * 1.380649e-23 * 298.15), 6)
1.0

entropy(T, N=6.02214076e+23)[source]#

Return the Sackur-Tetrode translational entropy.

\[S = Nk_B\left[\ln\left(\frac{q_{\text{trans}}}{N}\right) + \frac{5}{2}\right]\]

(O. Sackur, Ann. Phys. 36, 958 (1911); H. Tetrode, Ann. Phys. 38, 434 (1912); Atkins & de Paula, Physical Chemistry, 11th ed., eq. 13.24.) See also the standalone convenience function sackur_tetrode_entropy() for computing this directly from P and T without building a container volume by hand.

Parameters:: N (float)

heat_capacity_v(T, N=6.02214076e+23)[source]#

Return the constant-volume heat capacity Cv(T) contributed by this mode.

Parameters:

T (float)
N (float)

Returns:

float – Heat capacity, in J/K.

internal_energy(T, N=6.02214076e+23)[source]#

Return \(U(T)-U(0) = Nk_BT^2(\partial\ln q/\partial T)_V\).

Parameters:

T (float) – Absolute temperature, in K.
N (float) – Number of molecules/modes.

Returns:

float – Energy, in J (relative to this mode’s own ground state).

value(T)[source]#

Return the partition function q(T) for a single molecule/mode.

Parameters:: T (float) – Absolute temperature, in K.
Returns:: float

class chemistrykit.statmech.VibrationalPartitionFunctionHarmonic(frequency)[source]#

Bases: PartitionFunction

The vibrational partition function of a single harmonic-oscillator mode.

\[q_{\text{vib}} = \frac{1}{1-e^{-\Theta_{\text{vib}}/T}}, \qquad \Theta_{\text{vib}} = \frac{h\nu}{k_B}\]

with energies measured from the \(v=0\) ground vibrational level (Atkins & de Paula, Physical Chemistry, 11th ed., Ch. 13.2d). Unlike RotationalPartitionFunctionLinear’s classical limit, this is an exact closed form within the harmonic approximation, valid at every temperature – it is the harmonic-potential assumption itself (rather than a classical/high-T truncation) that is approximate, since real bonds are anharmonic (see chemistrykit.md.systems.pair_potentials.Morse). heat_capacity_v() is the Einstein-solid heat-capacity formula (A. Einstein, Ann. Phys. 22, 180 (1907)), reused here for a single vibrational mode rather than a 3D solid’s lattice vibrations.

Parameters:: frequency (float) – Vibrational frequency, in Hz.

Examples

The Einstein/vibrational heat capacity interpolates between 0 (T=0, the mode is frozen out) and the classical equipartition value \(Nk_B\) (\(T\to\infty\), one quadratic kinetic + one quadratic potential term):

>>> q = VibrationalPartitionFunctionHarmonic(frequency=8.7e13)  # ~HCl-like, ~2900 cm^-1
>>> round(float(q.heat_capacity_v(1.0, N=1.0)), 30)
0.0
>>> ratio = q.heat_capacity_v(1.0e6, N=1.0) / 1.380649e-23
>>> bool(0.99 < ratio < 1.0)
True

entropy(T, N=6.02214076e+23)[source]#

Return the entropy S(T) contributed by this mode.

Parameters:

T (float)
N (float)

Returns:

float – Entropy, in J/K.

classmethod from_wavenumber(wavenumber_cm_inv)[source]#

Build from a vibrational wavenumber in cm^-1 via \(\nu=c\tilde\nu\).

Parameters:: wavenumber_cm_inv (float) – Wavenumber \(\tilde\nu\), in cm^-1 (the usual IR-spectroscopy unit).
Return type:: VibrationalPartitionFunctionHarmonic
Returns:: VibrationalPartitionFunctionHarmonic

Examples

>>> q = VibrationalPartitionFunctionHarmonic.from_wavenumber(2886.0)  # HCl fundamental
>>> round(q.frequency / 1.0e13, 3)
8.652

heat_capacity_v(T, N=6.02214076e+23)[source]#

Return the constant-volume heat capacity Cv(T) contributed by this mode.

Parameters:

T (float)
N (float)

Returns:

float – Heat capacity, in J/K.

internal_energy(T, N=6.02214076e+23)[source]#

Return \(U(T)-U(0) = Nk_BT^2(\partial\ln q/\partial T)_V\).

Parameters:

T (float) – Absolute temperature, in K.
N (float) – Number of molecules/modes.

Returns:

float – Energy, in J (relative to this mode’s own ground state).

value(T)[source]#

Return the partition function q(T) for a single molecule/mode.

Parameters:: T (float) – Absolute temperature, in K.
Returns:: float

property vibrational_temperature: float#: The vibrational temperature \(\Theta_{\text{vib}}=h\nu/k_B\), in K.

chemistrykit.statmech.ising_2d_critical_temperature(coupling)[source]#

Return the exact square-lattice critical temperature \(T_c=2J/[k_B\ln(1+\sqrt2)]\).

Parameters:: coupling (float) – Exchange coupling J, in J.
Return type:: float
Returns:: float – Critical temperature, in K.

Examples

In units of \(J/k_B\), \(T_c\approx2.269\):

>>> round(ising_2d_critical_temperature(1.380649e-23), 4)
2.2692

chemistrykit.statmech.kramers_wannier_dual_coupling(K)[source]#

Return the Kramers-Wannier dual coupling \(K^*=-\tfrac12\ln\tanh K\).

The square-lattice Ising partition function at reduced coupling \(K=J/k_BT\) equals (up to an analytic prefactor) the one at \(K^*\), equivalently \(\sinh 2K\,\sinh 2K^*=1\): high temperature maps onto low temperature.

Parameters:: K (float or array_like) – Reduced coupling \(J/k_BT\) (positive).
Returns:: float or numpy.ndarray

Examples

The map is an involution, and its fixed point is the critical coupling \(K_c=\tfrac12\ln(1+\sqrt2)\):

>>> round(float(kramers_wannier_dual_coupling(kramers_wannier_dual_coupling(0.3))), 12)
0.3
>>> Kc = 0.5 * np.log(1.0 + np.sqrt(2.0))
>>> round(float(kramers_wannier_dual_coupling(Kc) - Kc), 12)
0.0

chemistrykit.statmech.ln_binomial(n, k)[source]#

Return \(\ln\binom{n}{k} = \ln n! - \ln k! - \ln(n-k)!\).

Parameters:

n (float or array-like of float)
k (float or array-like of float)

Returns:

float or ndarray

Examples

Agrees with direct combinatorics for small, exact integers:

>>> import numpy as np
>>> round(float(np.exp(ln_binomial(5, 2))), 6)
10.0

And, for large n, the peak degeneracy at \(k=n/2\) approaches Stirling’s classic large-n estimate \(\ln\binom{n}{n/2}\approx n\ln2\) (the basis of the large-M limit checked in chemistrykit.statmech.systems.lattice_gas):

>>> n = 100_000
>>> relative_error = abs(ln_binomial(n, n // 2) - n * np.log(2)) / (n * np.log(2))
>>> bool(relative_error < 1e-4)
True

chemistrykit.statmech.ln_factorial(n)[source]#

Return \(\ln(n!)\), via the exact log-gamma function \(\ln\Gamma(n+1)\).

Traditional statistical-mechanics derivations of combinatorial entropy approximate \(\ln n!\) with Stirling’s approximation (\(\ln n! \approx n\ln n - n\)); chemistrykit instead uses SciPy’s numerically exact gammaln (accurate to floating-point precision for every \(n\geq0\), including small n where Stirling’s approximation is poor), so no accuracy is sacrificed for the sake of matching the textbook derivation – the two agree closely only in the large-n limit that ln_binomial()’s docstring checks explicitly.

Parameters:: n (float or array-like of float) – Need not be an integer (\(\Gamma(n+1)\) generalizes n!).
Returns:: float or ndarray

Examples

>>> import numpy as np
>>> round(float(ln_factorial(5)), 6) == round(float(np.log(120)), 6)
True

chemistrykit.statmech.mean_speed(mass, temperature, k_b=1.380649e-23)[source]#

Mean speed, \(\bar v=\sqrt{8k_BT/(\pi m)}\).

Parameters:

mass (float)
temperature (float)
k_b (float, default chemistrykit.constants.K_B)

Returns:

float

chemistrykit.statmech.most_probable_speed(mass, temperature, k_b=1.380649e-23)[source]#

Most probable speed, \(v_p=\sqrt{2k_BT/m}\) – where \(f(v)\) peaks.

Parameters:

mass (float)
temperature (float)
k_b (float, default chemistrykit.constants.K_B)

Returns:

float

chemistrykit.statmech.rms_speed(mass, temperature, k_b=1.380649e-23)[source]#

Root-mean-square speed, \(v_{\text{rms}}=\sqrt{3k_BT/m}\).

Parameters:

mass (float)
temperature (float)
k_b (float, default chemistrykit.constants.K_B)

Returns:

float

chemistrykit.statmech.sackur_tetrode_entropy(mass, T, P, R_gas=8.31446261815324)[source]#

Molar translational entropy of an ideal gas via the Sackur-Tetrode equation.

A convenience wrapper around TranslationalPartitionFunction.entropy() that takes the pressure directly, using the ideal-gas volume per molecule \(V/N=k_BT/P\).

Parameters:

mass (float) – Molecular mass, in kg.
T (float) – Absolute temperature, in K.
P (float) – Pressure, in Pa.
R_gas (float)

Returns:

float – Molar entropy, in J/(mol K).

Examples

The textbook value for argon gas at 298.15 K and 1 bar is \(S_m^\circ\approx154.8\) J/(mol K) (Atkins & de Paula, Physical Chemistry, 11th ed., Table 13.1):

>>> import scipy.constants as sc
>>> S = sackur_tetrode_entropy(mass=39.948 * sc.atomic_mass, T=298.15, P=1.0e5)
>>> round(float(S), 1)
154.8

chemistrykit.statmech.second_virial_coefficient(pair_potential, T, length_scale, breakpoints=())[source]#

Return the molar second virial coefficient \(B_2(T)\) of a spherical pair potential.

Parameters:

pair_potential (Callable[[float], float]) – u(r) in J for a separation r in m; may return numpy.inf inside a hard core.
T (float) – Absolute temperature, in K.
length_scale (float) – A characteristic molecular size (e.g. \(\sigma\)), in m, used to split the integration range.
breakpoints (Sequence[float]) – Separations (m) where u is discontinuous (e.g. a square well’s edges), passed on to the quadrature.

Return type:

float

Returns:

float – \(B_2\), in m^3/mol (negative when attraction dominates).

Examples

Hard spheres give exactly \(\frac{2\pi}{3}N_A\sigma^3\):

>>> sigma = 3.4e-10
>>> hard_sphere = lambda r: np.inf if r < sigma else 0.0
>>> B2 = second_virial_coefficient(hard_sphere, 300.0, sigma, breakpoints=[sigma])
>>> round(B2 / (2 * np.pi / 3 * 6.02214076e23 * sigma**3), 8)
1.0

chemistrykit.statmech.thermal_de_broglie_wavelength(mass, temperature)[source]#

Return the thermal de Broglie wavelength \(\Lambda=h/\sqrt{2\pi mk_BT}\).

Parameters:

mass (float) – Particle mass, in kg.
temperature (float) – Absolute temperature, in K.

Returns:

float – Wavelength, in m.

Examples

For argon at room temperature, Lambda is a small fraction of an angstrom – far shorter than the mean interparticle spacing in a gas at atmospheric pressure, confirming that treating translational motion classically (as the ideal gas law does) is an excellent approximation there:

>>> import scipy.constants as sc
>>> wavelength = thermal_de_broglie_wavelength(mass=39.948 * sc.atomic_mass, temperature=298.15)
>>> bool(wavelength < 1e-10)
True