J. Chem. Theory Comput. 2013, 9, 3165−3169

我必须把对这篇文献的吐槽写在最前面

这学期学了统计热力学，前几天期末考试考完后就开始看这篇文献的。不得不说这篇文献在某些地方非常地坑爹，它的 Supporting Information 中有很多笔误。我在 follow 的时候必须充分翻阅各种资料才能确保公式的正确性。虽然后面的公式推导的思路是根据文献来的，但是我写的绝对比原本的要更清晰一些。

公式推导

自由电子气的熵

将电子作为量子 Fermi 气体，其巨正则配分函数为：

$$Z_G=\prod _p[1+{\mathrm{e}}^{-\beta ({\epsilon }_p-\mu )}]$$

式中 $\beta =1/kT$，${\epsilon }_p$ 代表第 $p$ 个可及能级，$\mu$ 为体系的化学势。由巨正则配分函数可得到巨热力学势

$$\mathrm{\Phi }=-\frac{1}{\beta }\ln Z_G=-\frac{1}{\beta }\sum _p\ln [1+{\mathrm{e}}^{-\beta ({\epsilon }_p-\mu )}]$$

将求和用积分近似（$J$ 是简并度，对电子而言等于 2）

$$\mathrm{\Phi }=-\frac{1}{\beta }\cdot (\frac{2\pi (2m{)}^{3/2}VJ}{h^3}{\int }_0^{\mathrm{\infty }}\ln [1+{\mathrm{e}}^{-\beta (\epsilon -\mu )}]{\epsilon }^{1/2}\mathrm{d}\epsilon )$$

定积分的部分可用分部积分得到

$$\begin{array}{rl}{\int }_0^{\mathrm{\infty }}\ln [1+{\mathrm{e}}^{-\beta (\epsilon -\mu )}]{\epsilon }^{1/2}\mathrm{d}\epsilon & =\frac{2}{3}(\ln [1+{\mathrm{e}}^{-\beta (\epsilon -\mu )}]\cdot {\epsilon }^{3/2}{|}_0^{\mathrm{\infty }}-{\int }_0^{\mathrm{\infty }}{\epsilon }^{3/2}\mathrm{d}(\ln [1+{\mathrm{e}}^{-\beta (\epsilon -\mu )}]))\\ & =\frac{2\beta }{3}{\int }_0^{\mathrm{\infty }}\frac{{\epsilon }^{3/2}}{1+{\mathrm{e}}^{\beta (\epsilon -\mu )}}\mathrm{d}\epsilon \end{array}$$

所以

$$\begin{array}{rl}\mathrm{\Phi }& =-\frac{2}{3}\frac{2\pi (2m{)}^{3/2}VJ}{h^3}{\int }_0^{\mathrm{\infty }}\frac{{\epsilon }^{3/2}}{1+{\mathrm{e}}^{\beta (\epsilon -\mu )}}\mathrm{d}\epsilon \\ & =-\frac{2}{3}\frac{JV{m}^{3/2}}{2^{1/2}{\pi }^2{\hslash }^3}{\int }_0^{\mathrm{\infty }}\frac{{\epsilon }^{3/2}}{1+{\mathrm{e}}^{\beta (\epsilon -\mu )}}\mathrm{d}\epsilon \end{array}$$

另一方面，已知

$$N=(\frac{\partial \mathrm{\Phi }}{\partial \mu })=\sum _{p}n({\epsilon }_p)$$

同样将求和用积分近似，有

$$\begin{array}{rl}N& =\sum _{p}n({\epsilon }_p)=\sum _p\frac{{\omega }_p}{1+{\mathrm{e}}^{\beta ({\epsilon }_p-\mu )}}\\ & =\frac{2\pi (2m{)}^{3/2}VJ}{h^3}{\int }_0^{\mathrm{\infty }}\frac{{\epsilon }^{1/2}}{1+{\mathrm{e}}^{\beta (\epsilon -\mu )}}\mathrm{d}\epsilon \\ & =\frac{JV{m}^{3/2}}{2^{1/2}{\pi }^2{\hslash }^3}{\int }_0^{\mathrm{\infty }}\frac{{\epsilon }^{1/2}}{1+{\mathrm{e}}^{\beta (\epsilon -\mu )}}\mathrm{d}\epsilon \end{array}$$

定义函数 $f_p$

$$f_p(a)=\frac{1}{p!}{\int }_0^{\mathrm{\infty }}\frac{x^p}{1+{\mathrm{e}}^x/a}\mathrm{d}x,a={\mathrm{e}}^{\mu /kT},x=\epsilon /kT$$

补充：

$$\begin{array}{rl}\mathrm{\Gamma }(\frac{1}{2}+n)& =(-\frac{1}{2}+n)!=\mathrm{\Pi }(-\frac{1}{2}+n)\\ & =\frac{(2n-1)!}{2^{2n-1}(n-1)!}\sqrt{\pi }\end{array}$$

$$(\frac{1}{2})!=\frac{\sqrt{\pi }}{2},(\frac{3}{2})!=\frac{3\sqrt{\pi }}{4}$$

所以摩尔巨热力学势为

$$\mathrm{\Phi }=-RT\frac{f_{3/2}(a)}{f_{1/2}(a)}$$

电子气的内能则为

$$U=\sum _p{\epsilon }_{p}n({\epsilon }_p)=-\frac{3}{2}\mathrm{\Phi }$$

熵由 Gibbs-Duhem 关系式得到

$$\begin{array}{rl}S& =\frac{U+pV}{T}-\frac{\mu N}{T}\\ & =\frac{5}{2}R\frac{f_{3/2}(a)}{f_{1/2}(a)}-R\ln a\end{array}$$

$$p=-(\frac{\partial \mathrm{\Phi }}{\partial V})=\frac{2}{3}\frac{J{m}^{3/2}}{2^{1/2}{\pi }^2{\hslash }^3}{\int }_0^{\mathrm{\infty }}\frac{{\epsilon }^{3/2}}{1+{\mathrm{e}}^{\beta (\epsilon -\mu )}}\mathrm{d}\epsilon$$

显然 $pV=-\mathrm{\Phi }=\frac{2}{3}U$。理想 Fermi 气体不遵循理想气体状态方程（$pV=NkT$）。

可推导得到温度为（默认取一个标准大气压）

$$T=[\frac{p^{\ominus }{h}^3}{Jk(2\pi mk{)}^{3/2}{f}_{3/2}(a)}{]}^{2/5}$$

强简并 Fermi 气体的情况（$T\ll T_F$）

Fermi 气体的能量为

$${\epsilon }_F=(\frac{6{\pi }^2}{J}{)}^{2/3}\frac{{\hslash }^2}{2m}{\rho }^{2/3}$$

$\rho =N/V$ 为数密度。以上公式的具体推导见 Wikipedia。

$$\mathrm{\Phi }(T)=-N{\epsilon }_F^{-3/2}{\int }_0^{\mathrm{\infty }}\frac{{\epsilon }^{3/2}}{1+{\mathrm{e}}^{\beta (\epsilon -\mu )}}\mathrm{d}\epsilon$$

$$1=\frac{3}{2}{\epsilon }_F^{-3/2}{\int }_0^{\mathrm{\infty }}\frac{{\epsilon }^{1/2}}{1+{\mathrm{e}}^{\beta (\epsilon -\mu )}}\mathrm{d}\epsilon$$

（中间涉及一些较复杂的过程，具体过程见文献或《统计热力学导论》P299-300）

$$\mu (T)={\epsilon }_F[1-\frac{{\pi }^2}{12}(\frac{kT}{{\epsilon }_F}{)}^2+O(T^4)]$$

$$\mathrm{\Phi }(T)=-\frac{2}{5}N{\epsilon }_F[1+\frac{5{\pi }^2}{12}(\frac{kT}{{\epsilon }_F}{)}^2+O(T^4)]$$

推导得到（摩尔）熵为

$$S=\frac{1}{T}(U+\frac{2}{3}U-\mu N)\approx \frac{{\pi }^2}{2}R\frac{T}{T_F}$$

其中定义了 Fermi 温度 $T_F={\epsilon }_F/k$。

热容

$$C_V=(\frac{\partial U}{\partial T})\approx \frac{{\pi }^2}{2}R\frac{T}{T_F}$$

$$C_p=T(\frac{\partial S}{\partial T})=\frac{{\pi }^2}{2}R\frac{T}{T_F}$$

完全非简并气体的情况（$T\gg T_F$）

这种情况下，$f_p(a)\approx a$。

$$S\approx \frac{5}{2}R-R\ln a$$

$$a\approx \frac{p^{\ominus }{h}^3}{Jk(2\pi mk{)}^{3/2}{T}^{5/2}}$$

得到

$$S=\frac{5}{2}R+R\ln [{\displaystyle \frac{J(2\pi mkT{)}^{3/2}kT}{p^{\ominus }{h}^3}}]$$

该表达式与使用电子的平动配分函数处理得到的结果是一致的。

迭代方法

文章中作者描述的方法是给定一个温度 $T_0$（即我们要得到该温度下电子的熵），以及 $a$ 的初猜。由该初猜可计算出相应的温度 $T$，将该温度与给定温度比较来判断 $a$ 的合理性。如果两者的差异大于要求的收敛限，则对 $a$ 进行校正，再一次计算温度。重复以上步骤直到温度 $T$ 收敛于要求的温度 $T_0$。这时利用 $a$ 便能算出熵。

但是作者在文中丝毫没有提到如何对 $a$ 进行校正，所以我使用了简单的二分法来求 $a$。另外 $f_p$ 函数计算中作者使用了 Romberg 方法进行数值积分，但是我并不会数值计算。从网上找的代码也难以完美用于这里的任务（被积分的函数有外部指定的参数），因此同样使用了简单的梯形公式来计算。所以整体的精度会差一些。求得不同温度下的熵后，可利用数值微分计算得到等压热容。并进一步算出焓等其它热力学量。同样地，我现在不会数值微分，因此其他的热力学量不方便求出。

结果

Fortran 代码上传到了GitHub。

用程序计算了 $251\mathrm{K}$ 到 $350\mathrm{K}$ 之间的电子的熵，部分结果如下（MathJax表格神烦，直接截图了）。

可以看出我得到的结果相较于文献计算的结果稍偏大。我估计自己得到的精度应该低于文献的，毕竟数值计算的部分用的是简单粗糙的方法。

虽然不能求得任意温度下的等压热容 $C_p$，但将熵对温度作图后发现，在至少 $251\sim 350\mathrm{K}$ 的范围内熵与温度近似为线性关系。进行线性拟合后得到斜率，利用公式 $$C_p=T(\frac{\partial S}{\partial T})$$ 可得到该温度范围内的等压热容。至于焓以及 Gibbs 自由能，则必须从 $0\mathrm{K}$ 开始计算等压热容的值，再根据 Kirchhoff 定律进行数值积分得到。

补充

把不同 $a$ 的情况下函数 $f_{3/2}$ 的图像画出来比较了一下

在 ${10}^{-1}$ 到 ${10}^4$ 这个数量级范围内，$f_p$在 $x=25$ 之前就基本趋于零了。所以我把积分的范围减小到了0至50（间隔数与之前一样为40000，之前的积分范围为0至120）。顺便把二分法的收敛限减小到 ${10}^{-9}$。重新计算后的结果有微小的改变，但是只在小数点第四位有变化。看来只有靠更好的数值方法才能有效地改善结果。

从一本Fortran科学计算的书上找了一个Simpson积分的代码拿来用，不过改进也不是特别大。另外，我想起来Origin可以做数值微分，所以算了 $273\sim 323\mathrm{K}$ 范围内的熵值。在Origin中算出等压热容 $C_p$，与文献值相比还是相对精确的。因为计算出的熵值与文献值相比是整体都偏大，但斜率则较接近。