.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/statmech/partition_functions/plot_03_gibbs_canonical_partition_function.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_statmech_partition_functions_plot_03_gibbs_canonical_partition_function.py: Gibbs's canonical ensemble: all thermodynamics from one partition function ============================================================================= Gibbs (1902) showed that for a system at fixed temperature every thermodynamic function follows from the canonical partition function :math:`Q=\sum_i e^{-E_i/k_BT}`: .. math:: A=-k_BT\ln Q,\qquad U=k_BT^2\frac{\partial\ln Q}{\partial T},\qquad S=\frac{U-A}{T},\qquad C_V=\frac{\partial U}{\partial T} Here the partition function of a vibrational mode is built by brute-force summation over its energy levels, and :math:`A`, :math:`U`, :math:`S`, :math:`C_V` are obtained from it by nothing but differentiation. They agree with the closed forms of :class:`~chemistrykit.statmech.VibrationalPartitionFunctionHarmonic`, and with :meth:`~chemistrykit.statmech.PartitionFunction.helmholtz_free_energy` (:math:`A=U-TS`). :class:`~chemistrykit.statmech.IdealGasMolecule` then combines independent modes, whose partition functions multiply. .. GENERATED FROM PYTHON SOURCE LINES 25-72 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from chemistrykit.constants import K_B from chemistrykit.statmech import IdealGasMolecule, VibrationalPartitionFunctionHarmonic mode = VibrationalPartitionFunctionHarmonic.from_wavenumber(1000.0) theta = mode.vibrational_temperature levels = K_B * theta * np.arange(400) # E_v = v h nu, measured from v = 0 def ln_Q(T): """Gibbs's sum over states, evaluated directly.""" return np.log(np.sum(np.exp(-levels / (K_B * T)))) T = np.linspace(100.0, 3000.0, 120) dT = 1e-2 * T lnQ = np.array([ln_Q(t) for t in T]) dlnQ = np.array([(ln_Q(t + d) - ln_Q(t - d)) / (2 * d) for t, d in zip(T, dT)]) A = -K_B * T * lnQ U = K_B * T**2 * dlnQ S = (U - A) / T Cv = np.gradient(U, T) fig, axes = plt.subplots(1, 2, figsize=(11, 4.5)) for values, closed, label, color in [ (A / (K_B * theta), mode.helmholtz_free_energy(T, N=1.0) / (K_B * theta), r"$A$", "crimson"), (U / (K_B * theta), mode.internal_energy(T, N=1.0) / (K_B * theta), r"$U$", "steelblue"), ]: axes[0].plot(T, closed, color=color, label=f"{label} closed form") axes[0].plot(T[::8], values[::8], "o", color=color, markersize=4, label=f"{label} from $\\ln Q$") axes[0].set_xlabel("T (K)") axes[0].set_ylabel(r"energy / $k_B\Theta_{\mathrm{vib}}$") axes[0].set_title("Free energy and internal energy from Q") axes[0].legend() axes[1].plot(T, mode.entropy(T, N=1.0) / K_B, color="seagreen", label=r"$S/k_B$ closed form") axes[1].plot(T[::8], S[::8] / K_B, "o", color="seagreen", markersize=4, label=r"$S = (U-A)/T$") axes[1].plot(T, mode.heat_capacity_v(T, N=1.0) / K_B, color="darkorange", label=r"$C_V/k_B$ closed form") axes[1].plot(T[::8], Cv[::8] / K_B, "o", color="darkorange", markersize=4, label=r"$C_V = \partial U/\partial T$") axes[1].set_xlabel("T (K)") axes[1].set_title("Entropy and heat capacity from Q") axes[1].legend() fig.tight_layout() .. image-sg:: /api/gallery/statmech/partition_functions/images/sphx_glr_plot_03_gibbs_canonical_partition_function_001.png :alt: Free energy and internal energy from Q, Entropy and heat capacity from Q :srcset: /api/gallery/statmech/partition_functions/images/sphx_glr_plot_03_gibbs_canonical_partition_function_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 73-76 The same machinery applies to a whole molecule: for independent modes the partition function factorizes, so ln q and every derived function are sums over modes. ``thermodynamic_functions`` bundles them: .. GENERATED FROM PYTHON SOURCE LINES 76-84 .. code-block:: Python co2_like = IdealGasMolecule(mass=7.3e-26, volume=0.0248, moment_of_inertia=7.2e-46, symmetry_number=2, vibrational_frequencies=[4.0e13, 2.0e13, 2.0e13, 7.0e13]) tf = co2_like.thermodynamic_functions(298.15) print(f"q = {tf.q:.3e}") print(f"U = {tf.U:.1f} J/mol, S = {tf.S:.2f} J/(mol K), Cv = {tf.Cv:.2f} J/(mol K), A = U - TS = {tf.A:.1f} J/mol") print(f"max |A_numeric - A_closed| / (k_B Theta) = {np.max(np.abs(A - mode.helmholtz_free_energy(T, N=1.0))) / (K_B * theta):.2e}") plt.show() .. rst-class:: sphx-glr-script-out .. code-block:: none q = 2.026e+33 U = 6888.0 J/mol, S = 213.81 J/(mol K), Cv = 28.83 J/(mol K), A = U - TS = -56858.0 J/mol max |A_numeric - A_closed| / (k_B Theta) = 9.09e-16 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.084 seconds) .. _sphx_glr_download_api_gallery_statmech_partition_functions_plot_03_gibbs_canonical_partition_function.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_03_gibbs_canonical_partition_function.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_03_gibbs_canonical_partition_function.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_03_gibbs_canonical_partition_function.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_