J. Chem. Theory Comput. 2013, 9, 3165−3169
我必须把对这篇文献的吐槽写在最前面
这学期学了统计热力学,前几天期末考试考完后就开始看这篇文献的。不得不说这篇文献在某些地方非常地坑爹,它的 Supporting Information 中有很多笔误。我在 follow 的时候必须充分翻阅各种资料才能确保公式的正确性。虽然后面的公式推导的思路是根据文献来的,但是我写的绝对比原本的要更清晰一些。
将电子作为量子 Fermi 气体,其巨正则配分函数为:
$$ Z_G=\prod\limits_p \Big[ 1+ {\mathrm e}^{-\beta(\varepsilon_p-\mu)} \Big] $$
式中 \(\beta=1/kT\),\(\varepsilon_p\) 代表第 \(p\) 个可及能级,\(\mu\) 为体系的化学势。由巨正则配分函数可得到巨热力学势
$$ \Phi = -\frac{1}{\beta}\ \ln{Z_G} = -\frac{1}{\beta} \sum\limits_p \ln{\Big[ 1+ {\mathrm e}^{-\beta(\varepsilon_p-\mu)} \Big] } $$
将求和用积分近似(\(J\) 是简并度,对电子而言等于 2)
$$ \Phi=-\frac{1}{\beta} \cdot \bigg( \frac{2\pi(2m)^{3/2}VJ}{h^3} \int_0^\infty \ln{\Big[ 1+ {\mathrm e}^{-\beta(\varepsilon-\mu)} \Big]} \varepsilon^{1/2} \, {\mathrm d}\varepsilon \bigg) $$
定积分的部分可用分部积分得到
$$ \begin{split} \int_0^\infty \ln{\Big[ 1+ {\mathrm e}^{-\beta(\varepsilon-\mu)} \Big]} \varepsilon^{1/2} \, {\mathrm d}\varepsilon &= \frac23 \bigg( \ln{\Big[ 1+ {\mathrm e}^{-\beta(\varepsilon-\mu)} \Big]}\cdot\varepsilon^{3/2} \Big|_0^\infty - \int_0^\infty \varepsilon^{3/2}\,{\mathrm d}\big( \ln{\Big[ 1+ {\mathrm e}^{-\beta(\varepsilon-\mu)} \big]} \Big) \bigg) \\ &= \frac{2\beta}{3} \int_0^\infty \frac{\varepsilon^{3/2}}{1+{\mathrm e}^{\beta(\varepsilon-\mu)}}\, {\mathrm d}\varepsilon \end{split} $$
所以
$$ \begin{split} {\Phi} &= -\frac23\ \frac{2\pi(2m)^{3/2}VJ}{h^3} \int_0^\infty \frac{\varepsilon^{3/2}}{1+{\mathrm e}^{\beta(\varepsilon-\mu)}}\, {\mathrm d}\varepsilon \\ &= -\frac23\ \frac{JVm^{3/2}}{2^{1/2}\pi^2\hbar^3} \int_0^\infty \frac{\varepsilon^{3/2}}{1+{\mathrm e}^{\beta(\varepsilon-\mu)}}\, {\mathrm d}\varepsilon \end{split} $$
另一方面,已知
$$ N = \bigg( \frac{\partial \Phi}{\partial \mu} \bigg) = \sum\limits_p n(\varepsilon_p) $$
同样将求和用积分近似,有
$$ \begin{split} N &= \sum\limits_p n(\varepsilon_p) = \sum\limits_p \frac{\omega_p}{1+{\mathrm e}^{\beta(\varepsilon_p-\mu)}} \\ &= \frac{2\pi(2m)^{3/2}VJ}{h^3} \int_0^\infty \frac{\varepsilon^{1/2}}{1+{\mathrm e}^{\beta(\varepsilon-\mu)}}\, {\mathrm d}\varepsilon \\ &= \frac{JVm^{3/2}}{2^{1/2}\pi^2\hbar^3} \int_0^\infty \frac{\varepsilon^{1/2}}{1+{\mathrm e}^{\beta(\varepsilon-\mu)}}\, {\mathrm d}\varepsilon \end{split} $$
定义函数 \(f_p\)
$$ f_p(a) = \frac{1}{p!} \int_0^\infty \frac{x^p}{1+{\mathrm e}^x/a}\,{\mathrm d}x ,\quad a={\mathrm e}^{\mu/kT} ,\quad x=\varepsilon/kT$$
补充:
$$ \begin{split} \Gamma \Big( \frac12+n \Big) &= \Big( -\frac12+n \Big)! = \Pi \Big( -\frac12+n \Big) \\ &= \frac{(2n-1)!}{2^{2n-1}(n-1)!}\sqrt{\pi} \end{split} $$
$$ \Big( \frac12 \Big)! = \frac{\sqrt{\pi}}{2},\quad \Big( \frac32 \Big)! = \frac{3\sqrt{\pi}}{4} $$
所以摩尔巨热力学势为
$$ \Phi = -RT\frac{f_{3/2}(a)}{f_{1/2}(a)} $$
电子气的内能则为
$$ U = \sum\limits_p \varepsilon_p n(\varepsilon_p) = -\frac32\ \Phi $$
熵由 Gibbs-Duhem 关系式得到
$$ \begin{split} {S} &= \frac{U+pV}{T} - \frac{\mu N}{T} \\ &= \frac52 R \frac{f_{3/2}(a)}{f_{1/2}(a)} - R \ln{a} \end{split} $$
$$ p = -\bigg( \frac{\partial \Phi}{\partial V} \bigg) = \frac23 \frac{Jm^{3/2}}{2^{1/2}\pi^2\hbar^3} \int_0^\infty \frac{\varepsilon^{3/2}}{1+{\mathrm e}^{\beta(\varepsilon-\mu)}}\, {\mathrm d}\varepsilon $$
显然 \( pV= -\Phi = \frac23 U\)。理想 Fermi 气体不遵循理想气体状态方程(\(pV=NkT\))。
可推导得到温度为(默认取一个标准大气压)
$$ T = \bigg[ \frac{p^\ominus h^3}{Jk(2\pi mk)^{3/2} f_{3/2}(a)} \bigg]^{2/5} $$
Fermi 气体的能量为
$$ \varepsilon_F = \bigg( \frac{6\pi^2}{J} \bigg)^{2/3} \frac{\hbar^2}{2m}\ \rho^{2/3} $$
\( \rho=N/V \) 为数密度。以上公式的具体推导见 Wikipedia。
$$ \Phi(T) = -N\varepsilon_F^{-3/2} \int_0^\infty \frac{\varepsilon^{3/2}}{1+{\mathrm e}^{\beta(\varepsilon-\mu)}}\, {\mathrm d}\varepsilon $$
$$ 1 = \frac32\ \varepsilon_F^{-3/2} \int_0^\infty \frac{\varepsilon^{1/2}}{1+{\mathrm e}^{\beta(\varepsilon-\mu)}}\, {\mathrm d}\varepsilon $$
(中间涉及一些较复杂的过程,具体过程见文献或《统计热力学导论》P299-300)
$$ \mu(T) = \varepsilon_F \bigg[ 1 - \frac{\pi^2}{12}\Big( \frac{kT}{\varepsilon_F} \Big)^2 + O(T^4) \bigg] $$
$$ \Phi(T) = -\frac25 N\varepsilon_F \bigg[ 1 + \frac{5\pi^2}{12}\Big( \frac{kT}{\varepsilon_F} \Big)^2 + O(T^4) \bigg] $$
推导得到(摩尔)熵为
$$ S = \frac{1}{T} \Big( U+\frac23 U -\mu N \Big) \approx \frac{\pi^2}{2} R \frac{T}{T_F} $$
其中定义了 Fermi 温度 \( T_F=\varepsilon_F/k \)。
热容
$$ C_V = \bigg( \frac{\partial U}{\partial T} \bigg) \approx \frac{\pi^2}{2} R \frac{T}{T_F} $$
$$ C_p = T \bigg( \frac{\partial S}{\partial T} \bigg) = \frac{\pi^2}{2} R \frac{T}{T_F} $$
这种情况下,\(f_p(a)\approx a\)。
$$ S \approx \frac52 R - R\ln{a} $$
$$ a \approx \frac{p^\ominus h^3}{Jk(2\pi mk)^{3/2} T^{5/2}} $$
得到
$$ S = \frac52 R + R\ln{ \bigg[ \dfrac{J(2\pi mkT)^{3/2} kT}{p^\ominus h^3} \bigg]} $$
该表达式与使用电子的平动配分函数处理得到的结果是一致的。
文章中作者描述的方法是给定一个温度 \(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 \bigg( \frac{\partial S}{\partial T} \bigg) $$ 可得到该温度范围内的等压热容。至于焓以及 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\),与文献值相比还是相对精确的。因为计算出的熵值与文献值相比是整体都偏大,但斜率则较接近。
[修改于 7年1个月前 - 2018/01/20 20:55:56]