求解电子的热力学量

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 \\ &\qquad\qquad – \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气体的情况(\(T\ll T_F\))

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} $$

完全非简并气体的情况(\(T\gg 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}\) 之间的电子的熵,部分结果如下。

\begin{array}{c|ccc}
\hline
T/\text{K} & S/{\mathrm J}\cdot\mathrm{K}^{-1}\cdot\mathrm{mol}^{-1} & S(\text{ref}) & A \\
\hline
260.00 & 20.4482 & 20.3635 & 1.711514 \\
261.00 & 20.5105 & 20.4256 & 1.691951 \\
262.00 & 20.5727 & 20.4876 & 1.672710 \\
263.00 & 20.6347 & 20.5495 & 1.653785 \\
264.00 & 20.6966 & 20.6113 & 1.635168 \\
265.00 & 20.7583 & 20.6729 & 1.616855 \\
266.00 & 20.8200 & 20.7344 & 1.598836 \\
267.00 & 20.8815 & 20.7957 & 1.581108 \\
268.00 & 20.9428 & 20.8569 & 1.563663 \\
269.00 & 21.0040 & 20.9180 & 1.546496 \\
270.00 & 21.0651 & 20.9790 & 1.529601 \\
271.00 & 21.1261 & 21.0398 & 1.512972 \\
272.00 & 21.1870 & 21.1005 & 1.496603 \\
273.00 & 21.2477 & 21.1611 & 1.480491 \\
274.00 & 21.3082 & 21.2216 & 1.464628 \\
275.00 & 21.3687 & 21.2819 & 1.449010 \\
276.00 & 21.4290 & 21.3421 & 1.433633 \\
277.00 & 21.4892 & 21.4022 & 1.418491 \\
278.00 & 21.5493 & 21.4621 & 1.403580 \\
279.00 & 21.6092 & 21.5219 & 1.388894 \\
280.00 & 21.6690 & 21.5816 & 1.374430 \\
\hline
\end{array}

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

虽然不能求得任意温度下的等压热容 \(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\),与文献值相比还是相对精确的。因为计算出的熵值与文献值相比是整体都偏大,但斜率则较接近。
\begin{array}{c|ccc}
\hline
T/\text{K} & C_p/{\mathrm J}\cdot\mathrm{K}^{-1}\cdot\mathrm{mol}^{-1} & C_p(\text{ref}) \\
\hline
273 & 16.5984 & 16.5202 \\
274 & 16.5496 & 16.5455 \\
275 & 16.6100 & 16.5707 \\
276 & 16.6704 & 16.5956 \\
277 & 16.6754 & 16.6204 \\
278 & 16.6800 & 16.6450 \\
279 & 16.6842 & 16.6694 \\
280 & 16.7440 & 16.6937 \\
\hline
\end{array}

最终版本

回家之后又研究了几天,做了如下改进。
把积分的部分换成了变步长Gauss Legendre积分,我对比了一下,和MatLab的数值积分的结果还是蛮接近的。

我觉得在积分这部分已经没法再奢求什么了。

优化求解 \(a\) 的部分曾尝试使用了Newton-Raphson法,但是效果并不好。比如迭代的过程中切线与x轴经常会相交于负半轴,这样就会导致后续的积分毫无意义(\(a\) 必须大于0)。所以依然沿用了二分法。

而为了保证二分法的效果,所以我参考了 J. Phys. Chem. 1994, 98, 6420-6424 中 \(a\) 的值用于辅助产生初猜。

上图中的曲线为拟合的 \(a\) 随温度 \(T\)变化的曲线,近似的规律是
$$ a(T) = 10^{A+\frac{B}{T+C}} $$
夹着曲线的两条直线为上下初猜的取值。

另外,我意识到之前统一把二分法的收敛限设为 \(10^{-9}\) 是不合理的。当 \(a\) 本身的数量级很大时,比如 \(10^{100}\),这时还要求 \(10^{-9}\) 程度的精确度是不对的。因为双精度的浮点数的有效位数是15位。所以现在对于温度小于 \(50\ {\mathrm K}\) 的情况,收敛限的数量级是动态变化的。变化规律是依据上图的。

现在的话,基本上对于温度大于 \(1\ {\mathrm K}\) 的情况,都能够算出来。

除此之外,还写了一个读取输入文件和写输出文件的功能,具体的使用方法写了一个Manual传到了GitHub。

发表评论

电子邮件地址不会被公开。 必填项已用*标注