Molecular Dynamics Simulation

From GPUMD
Jump to navigation Jump to search

告读者:本文还在编辑中。

  • 与本书配套的计算机程序将在此发布:https://github.com/brucefan1983/Molecular-Dynamics-Simulation
  • 2020 年 6 月 18 日更新:第 1 章(牛顿力学)已修订完毕,读者可开始阅读并找错了。本章无源代码。
  • 2020 年 6 月 18 日更新:第 2 章(分析力学)已修订完毕,读者可开始阅读并找错了。本章无源代码。
  • 2020 年 6 月 20 日更新:第 3 章(热力学)已修订完毕,读者可开始阅读并找错了。本章无源代码。
  • 2020 年 7 月 3 日更新:第 4 章(经典统计力学)已修订完毕,读者可开始阅读并找错了。本章无源代码。
  • 欢迎读者指出书稿中的谬误,提出修改意见:
    • 可加入 QQ 群 887975816 提出修改意见
    • 也可以给我发邮件:brucenju(at)gmail.com


  • 已收到的修改意见:
    • 张博 (华科)
      • [math]\vec{x}[/math] 改成 [math]\vec{r}[/math]?(等以后确定)
      • 热熔应该改成热容。(已改正)
      • 等体热容改成等容热容(已改正)
      • 热量和功的增量不要用全微分符号。(已改正)

Contents

牛顿力学 (已修订)

质点力学

运动学

  • 理论上,Newton力学里面最基本的研究对象是质点 (mass point),它是指一个有特定的质量 (mass) 但其大小对所研究的问题不重要的物体。有人也称它为一个粒子 (particle)。 一个粒子在三维空间的运动由一个位置函数 (position function) [math]\vec{x}(t)[/math] 完全描述。也就是说,在一个时刻 [math]t[/math],粒子的位置由一个有三个分量的矢量

$$ \vec{x}(t)=x(t)\vec{e}_x + y(t) \vec{e}_y + z(t) \vec{e}_z $$ 给出。当然,如果不指定一个参考系 (reference frame),我们将无法确定这三个分量的值。参考系的数学本质是坐标系 (coordinate system)。我们说粒子的位置函数完全地描述了粒子的运动性质,是因为其它的运动性质,比如速度和加速度,都可以由位置函数导出。

  • 速度 (velocity) 函数定义为位置函数对时间的一阶导数:

\begin{equation} \vec{v}(t) = \frac{d \vec{x}(t)}{d t} = \dot{\vec{x}}(t) =\frac{dx}{dt} \vec{e}_x + \frac{dy}{dt} \vec{e}_y + \frac{dz}{dt} \vec{e}_z. \end{equation}

  • 加速度 (acceleration) 函数 [math]\vec{a}(t)[/math] 为位置函数对时间的二阶导数,或者速度函数对时间的一阶导数:

\begin{equation} \vec{a}(t) = \frac{d ^2\vec{x}(t)}{d t^2} = \frac{d \vec{v}(t)}{d t} = \ddot{\vec{x}}(t) = \dot{\vec{v}}(t). \end{equation}

  • 以上对粒子运动的描述就是所谓的运动学

动力学

  • 动力学(dynamics)研究物体运动及其改变的原因。牛顿运动定律让我们可以在一定的条件下确定一个粒子的位置函数 [math]\vec{x}(t)[/math]
  • 牛顿第一定律。一个粒子在施加于它的合外力 [math]\vec{F}[/math] 为零时保持静止或者一个恒定的速度 [math]\vec{v}[/math]:

\begin{equation} \vec{F} = \vec{0} \Rightarrow \vec{v} = \text{const}. \end{equation} 该定律也叫惯性定律,或者惰性定律,意思是在合外力为零的情况下,一个粒子懒得改变它原来的运动状态:原来静止的话就继续静止;原来运动的话就继续原来的运动。

  • 牛顿第二定律。一个粒子的动量 (momentum) [math]\vec{p} [/math] 的时间变化率 [math]d\vec{p}/dt[/math] 正比于作用于它的合外力 [math]\vec{F}[/math]

\begin{equation} \vec{F} = \frac{d \vec{p}}{d t}. \end{equation} 一个粒子的动量定义为其惯性质量 [math]m[/math] 与速度 [math]\vec{v}[/math] 的乘积: $$ \vec{p}=m\vec{v} $$ 如果粒子的惯性质量是个常数的话,那么上式就变成 \begin{equation} \vec{F} = m\frac{d \vec{v}}{d t} = m \vec{a}. \end{equation} 我们强调了这里的 [math]m[/math] 是惯性质量,原因是牛顿力学里面还有另外一个质量,叫做引力质量。读者如果学过广义相对论的话,会知道惯性质量与引力质量是等价的。

  • 牛顿第三定律。每一个作用力都有一个与之大小相等、方向相反的反作用力。例如,如果有一个由粒子 2 作用在粒子 1 上的作用力 [math]\vec{F}_{12}[/math],那么一定同时存在一个由粒子 1 作用在粒子 2 上的反作用力 [math]\vec{F}_{21}[/math],它们大小相等,方向相反:

\begin{equation} \vec{F}_{12} = -\vec{F}_{21}. \end{equation} 当然,作用力和反作用力是相对的;把一对作用力中的哪一个叫做作用力,哪一个叫做反作用力是随意的。上式表达的是牛顿第三定律的弱形式 (weak form)。还有一个强形式 (strong form) 的牛顿第三定律: \begin{equation} \vec{F}_{12} \propto \vec{r}_{12} \equiv \vec{r}_{2} - \vec{r}_{1} \end{equation} 显然,满足了牛顿第三定律的强形式,就一定满足了牛顿第三定律的弱形式,但反之不然。

动量和动量守恒定律

  • 我们在陈述牛顿第二定律时已经定义了动量。根据牛顿第二定律,如果作用在一个粒子上的合力为零,那么该粒子的动量将保持为一个常数:

\begin{equation} \vec{F} = \vec{0} \Rightarrow \vec{p} = \text{常数}. \end{equation} 这就叫做动量守恒定律

角动量和角动量守恒定律

  • 角动量(angular momentum)被定义为:

\begin{equation} \vec{L} = \vec{x} \times \vec{p}. \end{equation} 因为位置 [math]\vec{x}[/math] 依赖于坐标系原点的选取,角动量也依赖于坐标原点的选择。类似地,我们定义力矩(torque or moment of force) [math]\vec{\tau}[/math]: \begin{equation} \vec{\tau} = \vec{x} \times \vec{F}. \end{equation} 其中,[math]\vec{F}[/math] 是作用在粒子上的力。让我们来考察角动量的时间变化率: \begin{equation} \frac{d \vec{L}}{dt} = \frac{d (\vec{x} \times \vec{p})}{dt} = \vec{v} \times \vec{p} + \vec{x} \times \vec{F} = \vec{x} \times \vec{F} = \vec{\tau}. \end{equation} 这就是说,角动量的变化率等于力矩。显然,如果作用在粒子上的力矩为零时,粒子的角动量是常数: \begin{equation} \vec{\tau} = \vec{0} \Rightarrow \vec{L} = \text{常数}. \end{equation} 这叫做角动量守恒定律。 细心的读者可能发现动量角动量两个词不大对称。确实,动量也可叫做线动量(linear momentum),与角动量相对,只不过我们通常将线动量简称为动量

动能和功

  • 除了动量,还可以定义一个仅仅与质量和速度有关的物理量,叫做动能(kinetic energy):

\begin{equation} T = \frac{1}{2}m \vec{v}^2 = \frac{1}{2}m \vec{v} \cdot \vec{v} = \frac{1}{2}m (v_x^2 + v_y^2 + v_z^2). \end{equation}

  • 动能如何得以改变呢?让我们来看看一个力对一个粒子做的(work)。在一个力 [math]\vec{F}[/math] 的作用下,如果粒子运动了一个微分位移 (differential displacement)

$$ d\vec{x} = dx \vec{e}_1 + dy \vec{e}_2 + dz \vec{e}_3 $$ 那么我们定义该力对该粒子做的微功 (differential work) 为: \begin{equation} d W = \vec{F} \cdot d\vec{x} = F_x dx + F_y dy + F_z dz. \end{equation}

  • 一般情况下,力是随时间变化的,可以写成时间的函数。如果在一个过程中,粒子受到一个力 [math]\vec{F}(t)[/math],而其位置的变化由轨道 [math]\vec{x}(t)~(t_1\lt t\lt t_2)[/math] 描述,那么我们说该过程中力 [math]\vec{F}[/math] 对粒子做的总功为如下积分:

\begin{equation} W = \int_{\vec{x}_1}^{\vec{x}_2} \vec{F} \cdot d\vec{x}. \end{equation} 其中,我们定义了 [math]\vec{x}_1 \equiv \vec{x}(t_1)[/math][math]\vec{x}_2 \equiv \vec{x}(t_2) [/math]。 力与微分位移的标量积可以变成如下形式: $$ \vec{F} \cdot d\vec{x} = \vec{F} \cdot \vec{v} dt = (\vec{F} dt) \cdot \vec{v} = d\vec{p} \cdot \vec{v} = m d\vec{v} \cdot \vec{v} = d \left( \frac{1}{2} m \vec{v}^2 \right). $$ 于是,我们有 \begin{equation} W = \frac{1}{2} m (\vec{v}_2^2 - \vec{v}_1^2) \equiv T_2 - T_1. \end{equation} 其中,我们定义了 [math]\vec{v}_i \equiv \vec{v}(t_i)[/math][math]T_i \equiv \frac{1}{2} m \vec{v}_i^2 ~(i=1, 2)[/math]

  • 以上推导告诉我们:在一个过程中,外力对一个粒子做的功等于粒子动能的改变量

势能

  • 在上面的讨论中,我们假定了一条连接起点和终点的具体的轨迹(即路径)。外力 [math]\vec{F}[/math] 作的功一般来说依赖于路径的选取。然而,如果外力做的功与具体的路径无关,而只与起点和终点有关,那么由矢量分析的定理可知,该外力可以写成一个标量场 [math]U(\vec{x})[/math] 的梯度的负值(这个负号是一种约定):

\begin{equation} \vec{F} = -\nabla U(\vec{x}). \end{equation} 该标量场 [math]U(\vec{x})[/math] 叫做粒子的势能场,简称为势能(potential energy),或者进一步简称为(potential)。与此对应的力称为保守力 (conservative force)。该保守力沿任意路径 [math]\vec{x}(t)~(t_1\lt t\lt t_2)[/math] 对粒子做的总功为: \begin{equation} W = - \int_{\vec{x}_1}^{\vec{x}_2} \nabla U(\vec{x}) \cdot d\vec{x} = - \int_{U_1}^{U_2} d U(\vec{x}) = U_1 - U_2. \end{equation} 其中,我们定义了 [math]U_1 = U(\vec{x}_1)[/math][math]U_2 = U(\vec{x}_2)[/math]。 将该式与前面得到的 [math]W = T_2 - T_1[/math] 比较可得: \begin{equation} T_1 + U_1 = T_2 + U_2. \end{equation} 这就是说,在一个保守力的作用下,任意过程中粒子的动能与势能的和都不改变。这个不改变的量称为机械能 (mechanical energy): \begin{equation} E = T + U. \end{equation} 所以,在保守力的作用下,粒子的机械能是守恒的。这就是机械能守恒原理。虽然称为原理,其实它并不是一个真正的原理,而是牛顿运动定律的推论。牛顿力学中真正的原理只有牛顿三大定律,其它的都是推论。

质点系力学

质点系的动量定理

  • 在具体讨论粒子系的动力学行为之前必须先弄清楚内力(internal force)和外力(external force)的区别。内力是系统中某个粒子作用于另一个粒子的,而外力来自于系统之外。内力满足牛顿第三定律。虽然外界施与粒子 [math]i[/math] 一个力 [math]\vec{F}^{\text{ext}}_{i}[/math] 的同时,粒子 [math]i[/math] 同时也施与外界一个大小相等、方向相反的反作用力,但由于我们的系统不包含外界,我们通常不会对外力利用第三定律。以上就是内力与外力的区别。
  • 对任意一个粒子 [math]i[/math],我们可以写下它的动力学方程:

\begin{equation} m_i \ddot{\vec{x}}_i = \sum_{j\neq i} \vec{F}_{ij} + \vec{F}^{\text{ext}}_{i}. \end{equation} 将上式左右两边都对指标 [math]i[/math] 求和,得 \begin{equation} \sum_i m_i \ddot{\vec{x}}_i = \sum_i \sum_{j\neq i} \vec{F}_{ij} + \sum_i \vec{F}^{\text{ext}}_{i}. \end{equation} 根据牛顿第三定律,等号右边的第一项等于零 (练习题)。系统的总质量定义为 [math]m = \sum_i m_i[/math]。如果定义一个“平均坐标” \begin{equation} \vec{x} = \frac{\sum_i m_i \vec{x}_i}{m}, \end{equation} 那么我们有 \begin{equation} m \ddot{\vec{x}} = \sum_i \vec{F}^{\text{ext}}_{i}. \end{equation} 这个式子看上去很像一个质量为系统总质量 [math]m[/math],坐标为“平均坐标” [math]\vec{x}[/math] 的“粒子”的动力学方程;该“粒子”所受的合外力为整个粒子系所受的合外力。我们称这个等效的“粒子”为粒子系的质心(center of mass)。质心的质量就是整个系统的质量;质心坐标就是上述“平均坐标”;质心的运动满足上述等效的牛顿第二定律。

  • 由质心坐标可以定义质心速度

\begin{equation} \dot{\vec{x}} = \frac{\sum_i m_i \dot{\vec{x}}_i}{m} \end{equation} 和质心动量 \begin{equation} \vec{p} = m \dot{\vec{x}} = \sum_i m_i \dot{\vec{x}}_i. \end{equation} 于是,质心的牛顿第二定律可用质心动量表达为: \begin{equation} \frac{d \vec{p}}{dt} = \sum_i \vec{F}^{\text{ext}}_{i}. \end{equation} 如果系统受到的合外力为零,那么系统的质心动量(即系统的总动量)是守恒的。这就是质点系的动量守恒定律。

质点系的角动量守恒定理

类似地,由牛顿第三定律可以证明内力对系统的总力矩的贡献也是零,即系统所受的总力矩等于外力的总力矩 [math]\vec{\tau}^{\text{ext}}[/math]: \begin{equation} \vec{\tau}^{\text{ext}} = \sum_i^N \vec{x}_i \times \vec{F}_i^{\text{ext}}. \end{equation} 如果定义系统的总角动量为 \begin{equation} \vec{L} = \sum_i^N \vec{L}_i = \sum_i^N \vec{x}_i \times \vec{p}_i, \end{equation} 那么可以证明如下的角动量定理: \begin{equation} \frac{d \vec{L}}{dt} = \vec{\tau}^{\text{ext}}. \end{equation} 这就是说,外力产生的总力矩等于系统总角动量的时间变化率。如果外力产生的总力矩等于零,则有系统的总角动量守恒。这就是质点系的角动量守恒定理。

质点系的功与能

不同于力和力矩,内力对系统中的粒子做的总功并不一定为零。例如,两个粒子在相互的排斥力的作用下从静止开始做反向的加速运动,它们之间的内力做了正功。在一个微小过程中,如果第 [math]i[/math] 个粒子的微分位移为 [math]d\vec{x}_i[/math],那么内力对系统中的粒子做的功为 \begin{equation} dW^{\text{int}} = \sum_i^N \vec{F}_i^{\text{int}} \cdot d\vec{x}_i, \end{equation} 外力对系统中的粒子做的功为 \begin{equation} dW^{\text{ext}} = \sum_i^N \vec{F}_i^{\text{ext}} \cdot d\vec{x}_i. \end{equation} 系统中的粒子得到的总功为 \begin{equation} dW = dW^{\text{ext}} + dW^{\text{int}}. \end{equation} 因为对每一个粒子,有 \begin{equation} (\vec{F}_i^{\text{int}} + \vec{F}_i^{\text{ext}}) \cdot d\vec{x}_i = d T_i, \end{equation} 所以,综合起来有 \begin{equation} dW=dT. \end{equation} 其中,[math]T[/math] 定义为系统的总动能: \begin{equation} T = \sum_i^N T_i = \sum_i^N \left( \frac{1}{2} m_i \vec{v}_i^2 \right). \end{equation} 这就是说,内力和外力做的功的总和等于系统总动能的改变。只有内力和外力做的总功为零时,系统的动能才守恒。


对于有相互作用的多粒子系统,如果内力和外力都是保守力,那么系统的总势能可以表达为 \begin{equation} U = \sum_i U_i + \frac{1}{2} \sum_i^N \sum_{j\neq i}^N U_{ij}, \end{equation} 其中,[math]U_i[/math] 是第 [math]i[/math] 个粒子在外力场中的势能,[math]U_{ij}[/math] 是系统中由 [math]i[/math][math]j[/math] 的相互作用导致的势能。

习题

简谐振子

  • 题目:

考虑一维系统。假设粒子在如下势能场中运动: \begin{equation} U(x) = \frac{1}{2} k x^2. \end{equation} 这里的 [math]k[/math] 为常数。

    • 请写出该粒子的运动方程。
    • 请自行假设一定的初始条件求解其运动轨迹。
    • 证明其运动过程中机械能守恒。
  • 解答:
    • 粒子在该势场中受到的力由势场的负梯度给出:

\begin{equation} F = - \frac{dU(x)}{dx} = -kx. \end{equation} 该力总是指向原点,且大小正比于偏离原点的距离,它表达了胡克定律。 于是,粒子的运动方程可写为: \begin{equation} m\ddot{x} = - kx. \end{equation}

    • 定义一个具有频率的量纲的量

\begin{equation} \omega = \sqrt{k/m}, \end{equation} 可将运动方程写为 \begin{equation} \ddot{x} + \omega^2 x = 0. \end{equation} 该方程的通解为 \begin{equation} x = B \cos(\omega t) + C \sin(\omega t) \end{equation} 两个待定系数 [math]B[/math][math]C[/math] 可以由两个初始条件确定: \begin{equation} x(t=0) = x_0 = B; \quad \dot{x}(t=0) = v_0 = \omega C. \end{equation} 所以,粒子的坐标函数为: \begin{equation} x = x_0 \cos(\omega t) + \frac{v_0}{\omega} \sin(\omega t) = A \sin(\omega t + \delta). \end{equation} 其中, $$A = \left(x_0^2 + v_0^2/\omega^2\right)^{1/2}$$ $$\tan \delta = B/C$$

    • 粒子的动能 [math]T[/math] 和势能 [math]U[/math]

\begin{equation} T = \frac{1}{2} m \dot{x}^2 = \frac{1}{2} m A^2 \omega^2 \cos^2(\omega t + \delta) = \frac{1}{2} k A^2 \cos^2(\omega t + \delta) , \end{equation} \begin{equation} U = \frac{1}{2} k x^2 = \frac{1}{2} k A^2 \sin^2(\omega t + \delta). \end{equation} 正如我们所期望的,该保守力系统的总能量 [math]E[/math] 是守恒的: \begin{equation} E = T + U = \frac{1}{2} k A^2. \end{equation}

质点系内力的和是零

  • 证明这个论断。

质点系内力对系统的总力矩的贡献是零

  • 证明这个论断。


质点系的角动量定理

  • 证明质点系的角动量定理:

\begin{equation} \frac{d \vec{L}}{dt} = \vec{\tau}^{\text{ext}}. \end{equation}

柯尼希 (König) 定理

  • 题目:

证明质点系的总动能可以写成质心动能 [math]T_c[/math] 和相对动能 [math]T_r[/math] 两部分的和: \begin{equation} T = \frac{1}{2} m \vec{v}^2 + \frac{1}{2}\sum_i^N m_i (\vec{v}_i-\vec{v})^2 \equiv T_c + T_r. \end{equation} 这个分解叫做柯尼希定理

  • 证明:

[math]\vec{v}_i[/math] 写成 [math]\vec{v} + (\vec{v}_i-\vec{v})[/math],我们有 \begin{equation} T = \frac{1}{2}\sum_i^N m_i \vec{v}_i^2 = \frac{1}{2}\sum_i^N m_i [(\vec{v} + (\vec{v}_i-\vec{v})]^2. \end{equation} 展开得 \begin{equation} T = \frac{1}{2}\sum_i^N m_i \vec{v}^2 + \frac{1}{2}\sum_i^N m_i (\vec{v}_i-\vec{v})^2 + \sum_i^N m_i \vec{v} \cdot (\vec{v}_i-\vec{v}). \end{equation} 可见,只要证明上式中的最后一项等于零就可以了。证明如下: \begin{equation} \sum_i^N m_i \vec{v} \cdot (\vec{v}_i-\vec{v}) =\sum_i^N m_i \vec{v} \cdot \vec{v}_i - \sum_i^N m_i \vec{v} \cdot \vec{v} = \vec{v} \cdot \vec{p} - m \vec{v} \cdot \vec{v} = 0. \end{equation}

分析力学 (已修订)

广义坐标

  • 我们知道,一个由 [math]N[/math] 个质点构成的力学系统需要用 [math]3N[/math] 个坐标分量及其相应的速度分量来描述其运动状态。我们说,该系统的力学自由度 (degree of freedom) 为 [math]3N[/math]。 但如果是一个由这 [math]N[/math] 个质点构成的刚体,则自由度只有 6:3 个平动自由度和 3 个转动自由度。刚体的自由度之所以小于自由质点系统的自由度,是因为刚体系统中有约束 (constraint)。
  • 约束的分类很复杂,但我们只关注最简单的一种:稳定的几何约束。考虑一个由 [math]N[/math] 个质点构成的力学系统,其中第 [math]i[/math] 个质点的位置为 [math]\vec{x}_i[/math]。一个稳定的几何约束一般可以表达为:

\begin{equation} f(\vec{x}_1, \vec{x}_2,\cdots) = 0. \end{equation} 其中,[math]f[/math] 是个一般的函数。

  • 容易看出,一个原本有 [math]M[/math] 个自由度的体系,若受到 [math]m[/math] 个几何约束,则其自由度为 [math]s=M-m[/math]。要描述该体系,我们不一定需要使用原来的 [math]M[/math] 个坐标和速度,而可以构造 [math]s[/math] 个新的坐标和速度。这 [math]s[/math] 个新的坐标不一定限于在原来的 [math]M[/math] 个老坐标中挑选,而完全可以另行选取。这样选取的新坐标叫做广义坐标 (generalized coordinate)。广义坐标不一定能三个一组地构成矢量,而且某个广义坐标也不一定具有长度的量纲。相应地,某个广义速度 (generalized velocity) 也不一定具有速度的量纲。

虚功原理和达朗伯 (d'Alembert) 原理

虚功原理

  • 在有约束的情况下,力学系统各个部分虽然不能随心所欲地运动,但总能作某种运动。我们把考虑约束之前的某个自由度在有约束的情况下所具有的一切可能的位移叫做虚位移 (virtual displacement),记为 [math]\delta \vec{x}_i[/math]
  • 系统中某个质点 [math]i[/math] 除了可能受到普通的力 [math]\vec{F}_i[/math] (叫做主动力) 之外,还可能受到约束反力 [math]\vec{N}_i[/math]。若该质点保持力学平衡,则有

\begin{equation} \vec{F}_i + \vec{N}_i = 0. \end{equation}

  • 正如力与位移的标量积定义为功,力与虚位移的标量积定义为虚功 (virtual work)。显然,质点 [math]i[/math] 平衡时主动力与约束反力的虚功之和 [math]\delta W_i[/math] 为零:

\begin{equation} \delta W_i = (\vec{F}_i + \vec{N}_i) \cdot \delta \vec{x}_i = 0. \end{equation} 对所有质点进行求和,则有 \begin{equation} \delta W = \sum_{i=1}^{N} \delta W_i = \sum_{i=1}^{N} (\vec{F}_i + \vec{N}_i) \cdot \delta \vec{x}_i = 0. \end{equation} 在很多重要的情况下,作为内力的各个约束反力做的总虚功为零: \begin{equation} \sum_{i=1}^{N} \vec{N}_i \cdot\delta \vec{x}_i = 0. \end{equation} 这种情况下的约束叫做理想约束。在理想约束作用下,我们有 \begin{equation} \sum_{i=1}^{N} \vec{F}_i \cdot \delta \vec{x}_i = 0. \end{equation} 该式表达的就是虚功原理 (principle of virtual work)。

达朗伯原理

以上考虑的是静力学问题。在动力学问题中,质量为 [math]m_i[/math]、坐标为 [math]\vec{x}_i[/math] 的质点的运动方程可以写为: \begin{equation} \vec{F}_i + \vec{N}_i = m_i\ddot{\vec{x}}_i. \end{equation} 若选取跟随质点一同运动的参考系,则上式可改写为 \begin{equation} \vec{F}_i - m_i\ddot{\vec{x}}_i + \vec{N}_i = 0. \end{equation} 其中,[math]- m_i\ddot{\vec{x}}_i[/math] 可以理解为惯性力。类似于对虚功原理的推导,我们有如下等式: \begin{equation} \sum_{i=1}^{N} (\vec{F}_i - m_i\ddot{\vec{x}}_i+ \vec{N}_i) \cdot \delta \vec{x}_i = 0. \end{equation} 在理想约束的作用下,我们有 \begin{equation} \sum_{i=1}^{N} \vec{F}_i \cdot \delta \vec{x}_i - \sum_{i=1}^{N} m_i\ddot{\vec{x}}_i \cdot \delta \vec{x}_i = 0. \end{equation} 该式表达的就是达朗伯原理。它是下一小节讨论的出发点。

拉格朗日 (Lagrange) 方程

  • 考虑 [math]N[/math] 个质点组成的力学系统。质点 [math]i[/math] 的质量和坐标分别为 [math]m_i[/math][math] \vec{x}_i[/math] 。其中 [math]i=1, 2, \cdots, N[/math]。 若有 [math]m[/math] 个理想的几何约束,则可以定义 [math]s=3N-m[/math] 个广义坐标 [math]q_1, q_2, \cdots q_s[/math]。通常将这些坐标的集合简记为 [math]q[/math]。原来的坐标可以表达为广义坐标的函数:

\begin{equation} \vec{x}_i = \vec{x}_i(q) \quad (i = 1, 2, \cdots, N). \end{equation}

  • 运用上述坐标变换,可以将达朗伯原理的表达式写成

\begin{equation} \vec{F}_i \cdot \frac{\partial \vec{x}_i}{\partial q_{\alpha}} \delta q_{\alpha} - m_i \ddot{\vec{x}_i} \cdot \frac{\partial \vec{x}_i} {\partial q_{\alpha}} \delta q_{\alpha} = 0. \end{equation} 注意:

    • 我们用了爱因斯坦求和约定。
    • 虚位移不是原来的坐标所特有的,新的广义坐标也可以有虚位移,因为一个坐标如果没有虚位移那就没有让它成为坐标的必要了。
  • 因为新的广义坐标是相互独立(因为广义坐标个数等于自由度个数,我们总可以做到这一点)的,所以根据上式我们有下式:

\begin{equation} \vec{F}_i \cdot \frac{\partial \vec{x}_i}{\partial q_{\alpha}} - m_i \ddot{\vec{x}_i} \cdot \frac{\partial \vec{x}_i} {\partial q_{\alpha}} = 0. \quad (\alpha = 1, 2, \cdots, s). \end{equation} 下面的任务就是从该式出发推导拉格朗日方程。

  • 上式中的加速度 [math]\ddot{\vec{x}_i}[/math] 是比较讨厌的,我们可以用求导的乘积率换掉它:

\begin{equation} m_i \ddot{\vec{x}_i} \cdot \frac{\partial \vec{x}_i} {\partial q_{\alpha}} = m_i \frac{d}{dt} \left( \dot{\vec{x}_i} \cdot \frac{\partial \vec{x}_i} {\partial q_{\alpha}} \right) -m_i \dot{\vec{x}_i} \cdot \frac{d}{dt} \left( \frac{\partial \vec{x}_i} {\partial q_{\alpha}} \right) \end{equation}


  • 上面的等式中出现了速度 [math]\dot{\vec{x}_i}[/math] 和位置 [math]\vec{x}_i[/math] 的混合,还是比较讨厌。你是不是想:要是

\begin{equation} \frac{\partial \vec{x}_i} {\partial q_{\alpha}} = \frac{\partial \dot{\vec{x}}_i}{\partial \dot{q}_{\alpha}}, \end{equation} \begin{equation} \frac{d}{dt} \left( \frac{\partial \vec{x}_i} {\partial q_{\alpha}} \right) = \frac{\partial \dot{\vec{x}}_i}{\partial q_{\alpha}} \end{equation} 就好了?这两个式子确实是成立的。对于第一个,我们证明右边等于左边: \begin{equation} \frac{\partial \dot{\vec{x}}_i}{\partial \dot{q}_{\alpha}} = \frac{\partial}{\partial \dot{q}_{\alpha}} \left( \frac{\partial \vec{x}_i} {\partial q_{\beta}} \dot{q}_{\beta}\right) = \frac{\partial \vec{x}_i} {\partial q_{\beta}} \delta_{\alpha\beta} = \frac{\partial \vec{x}_i} {\partial q_{\alpha}}. \end{equation} 对于第二个,我们证明左边等于右边: \begin{equation} \frac{d}{dt} \left( \frac{\partial \vec{x}_i} {\partial q_{\alpha}} \right) = \frac{\partial}{\partial q_{\beta}} \left( \frac{\partial \vec{x}_i} {\partial q_{\alpha}} \right) \dot{q}_{\beta} = \frac{\partial}{\partial q_{\alpha}} \left( \frac{\partial \vec{x}_i} {\partial q_{\beta}} \right) \dot{q}_{\beta} = \frac{\partial \dot{\vec{x}}_i}{\partial q_{\alpha}} \end{equation} 利用它们,我们马上有: \begin{align} m_i \ddot{\vec{x}_i} \cdot \frac{\partial \vec{x}_i} {\partial q_{\alpha}} =& m_i \frac{d}{dt} \left( \dot{\vec{x}_i} \cdot \frac{\partial \dot{\vec{x}}_i}{\partial \dot{q}_{\alpha}} \right) -m_i \dot{\vec{x}_i} \cdot \frac{\partial \dot{\vec{x}}_i}{\partial q_{\alpha}} \nonumber \\ =& \frac{d}{dt} \left[ \frac{\partial}{\partial \dot{q}_{\alpha}} \left( \frac{m_i\dot{\vec{x}_i} \cdot \dot{\vec{x}}_i}{2}\right) \right] -\frac{\partial}{\partial q_{\alpha}} \left( \frac{m_i\dot{\vec{x}_i} \cdot \dot{\vec{x}}_i}{2}\right) \end{align}


  • 神奇的时刻到了:上式中第二个等号右边的两对小括号中的表达式都是系统的动能 [math]T[/math]。于是,我们最终可以将达朗伯原理的表达式写成:

\begin{equation} \frac{d}{dt} \left(\frac{\partial T}{\partial \dot{q}_{\alpha}}\right) -\frac{\partial T}{\partial q_{\alpha}} =\vec{F}_i \cdot \frac{\partial \vec{x}_i}{\partial q_{\alpha}}. \quad (\alpha = 1, 2, \cdots, s) \end{equation} 这就是我们要推导的拉格朗日方程。右边的表达式叫做广义力 (generalized force)。对于保守力,主动力的虚功可以写成 \begin{equation} \vec{F}_i \cdot \frac{\partial \vec{x}_i} {\partial q_{\alpha}} \delta q_{\alpha} = \vec{F}_i \cdot \delta \vec{x}_i = -\delta V. \end{equation} 其中,[math]V[/math] 是体系的势能函数。于是, \begin{equation} \vec{F}_i \cdot \frac{\partial \vec{x}_i} {\partial q_{\alpha}} = - \frac{ \delta V} { \delta q_{\alpha} } = - \frac{ \partial V} { \partial q_{\alpha} }. \end{equation} 所以,对于保守力体系,拉格朗日方程可写为: \begin{equation} \frac{d}{dt} \left(\frac{\partial T}{\partial \dot{q}_{\alpha}}\right) -\frac{\partial T}{\partial q_{\alpha}} = - \frac{ \partial V} { \partial q_{\alpha} }. \quad (\alpha = 1, 2, \cdots, s) \end{equation} 势能 [math]V[/math] 仅仅是广义坐标的函数,而不是广义速度 [math]\dot{q}_{\alpha}[/math] 的函数, 故有 \begin{equation} \frac{d}{dt} \left(\frac{\partial (T-V)}{\partial \dot{q}_{\alpha}}\right) -\frac{\partial (T-V)}{\partial q_{\alpha}} = 0. \quad (\alpha = 1, 2, \cdots, s) \end{equation} 这里出现的动能与势能的差是一个很重要的量,叫做拉格朗日量 (Lagrangin),记为 \begin{equation} L(q, \dot{q}) = T(\dot{q}) - V(q). \end{equation} 用拉格朗日量可以将拉格朗日方程写成更加简洁的形式: \begin{equation} \boxed{ \frac{d}{dt} \left(\frac{\partial L}{\partial \dot{q}_{\alpha}}\right) -\frac{\partial L}{\partial q_{\alpha}} = 0. \quad (\alpha = 1, 2, \cdots, s) } \end{equation}

哈密顿 (Hamilton) 方程

  • 由于拉格朗日量具有能量的量纲,故当广义坐标具有长度的量纲时,偏导数 [math]\frac{\partial L}{\partial \dot{q}_{\alpha}}[/math] 具有动量的量纲。我们称这个偏导数为广义动量 (generalized momentum),记为

\begin{equation} \boxed{p_{\alpha} = \frac{\partial L}{\partial \dot{q}_{\alpha}}}. \end{equation} 这个定义是要牢记的。注意,当广义坐标不具有长度的量纲时,相应的广义动量也不具有动量的量纲。运用广义动量可以将拉格朗日方程写成更加简洁的形式: \begin{equation} \boxed{ \dot{p}_{\alpha} -\frac{\partial L}{\partial q_{\alpha}} = 0. \quad (\alpha = 1, 2, \cdots, s) } \end{equation} 如果你也能记住这个公式那就太好了。

  • 虽然广义动量出现在拉格朗日方程中,但我们要清楚的是拉格朗日量是广义坐标和广义速度的函数,而不是广义动量的函数。用勒让德变换可以改变一个多元函数的独立变量。为此,我们先写下拉格朗日量的全微分:

\begin{equation} d L(q, \dot{q}) = \frac{\partial L}{\partial q_{\alpha}} d q_{\alpha} + \frac{\partial L}{\partial \dot{q}_{\alpha}} d \dot{q}_{\alpha} = \dot{p}_{\alpha} d q_{\alpha} + p_{\alpha} d \dot{q}_{\alpha}. \end{equation} 下面的拉格朗日变换能将独立变量从广义速度变成广义动量: \begin{equation} L \rightarrow p_{\alpha} \dot{q}_{\alpha} - L. \end{equation} 该式箭头右边的表达式是一个很重要的量,记为 \begin{equation} H= p_{\alpha} \dot{q}_{\alpha} - L. \end{equation} 该量称为哈密顿量 (Hamiltonian)。我们来看看它是不是真的是广义动量的函数了。为此,我们计算它的全微分: \begin{equation} d H = p_{\alpha} d\dot{q}_{\alpha} + \dot{q}_{\alpha} d p_{\alpha} - (\dot{p}_{\alpha} d q_{\alpha} + p_{\alpha} d \dot{q}_{\alpha}) =\dot{q}_{\alpha} d p_{\alpha} - \dot{p}_{\alpha} d q_{\alpha}. \end{equation} 看来哈密顿量确实是广义动量以及广义坐标的函数,而不是广义速度的函数了。

  • 既然哈密顿量是广义动量和广义坐标的函数,那么根据全微分的定义,我们又有

\begin{equation} d H = \frac{\partial H}{\partial q_{\alpha}} d q_{\alpha} + \frac{\partial H}{\partial p_{\alpha}} d p_{\alpha}. \end{equation} 对比以上两个式子,我们得到如下两组重要的方程: \begin{equation} \boxed{\dot{q}_{\alpha} = \frac{\partial H}{\partial p_{\alpha}}, \quad (\alpha = 1 , 2, \cdots, s)} \end{equation} \begin{equation} \boxed{\dot{p}_{\alpha} = - \frac{\partial H}{\partial q_{\alpha}}. \quad (\alpha = 1 , 2, \cdots, s)} \end{equation} 这 [math]2s[/math] 个一阶微分方程组称为哈密顿正则方程 (canonical equations)。这里的“正则”是“简单而且对称”的意思,不要被这个怪怪的形容词吓到了。我相信你能够记住这两组方程,因为它们太重要了。

相空间

广义坐标和广义动量是相互独立的变量。我们可以将 [math]s[/math] 个广义坐标 [math]\{q_{\alpha}\}_{\alpha=1}^{s}[/math][math]s[/math] 个广义动量 [math]\{p_{\alpha}\}_{\alpha=1}^{s}[/math] 看成一个 [math]2s[/math] 维“空间”的“坐标”。这个抽象的“空间”叫做相空间 (phase space)。一组给定的广义坐标和广义动量叫做相空间的一个相点 (phase point)。另外,根据哈密顿正则方程,只要给定一个初始条件,即初始时刻的各个广义坐标和广义动量,就可以唯一地确定任意时刻的各个广义坐标和广义动量,即如下 [math]2s[/math] 个函数: \begin{equation} q_{\alpha} = q_{\alpha}(t), \quad p_{\alpha} = p_{\alpha}(t). \quad (\alpha = 1, 2, \cdots, s) \end{equation} 这 [math]2s[/math] 个函数将在 [math]2s[/math] 维的相空间给出一条轨迹,叫做相轨迹 (phase trajectory)。随着时间的推移,一条相轨迹会在相空间跑动,历经很多不同的相点。一个系统中不同的初始条件会给出不同的相轨迹,而两条不同的相轨迹绝对不会相交于某一个相点。

泊松 (Poisson) 括号和刘维尔算符

一个一般的物理量可以表达为相空间坐标的函数 $$ A=A(q_{\alpha}, p_{\alpha}) $$ 我们求它的时间导数 $$ \frac{dA}{dt}=\frac{\partial A}{\partial q_{\alpha} } \dot{q}_{\alpha} +\frac{\partial A}{\partial p_{\alpha} } \dot{p}_{\alpha} $$ 利用哈密顿正则方程,可得 $$ \frac{dA}{dt}=\frac{\partial A}{\partial q_{\alpha} } \frac{\partial H}{\partial p_{\alpha} } -\frac{\partial A}{\partial p_{\alpha} } \frac{\partial H}{\partial q_{\alpha} } $$


定义任意两个物理量之间的泊松括号 $$ \{A,B\}=\frac{\partial A}{\partial q_{\alpha} } \frac{\partial B}{\partial p_{\alpha} } -\frac{\partial A}{\partial p_{\alpha} } \frac{\partial B}{\partial q_{\alpha} } $$ 我们有 $$ \frac{dA}{dt}= \{A,H\} $$


一个物理量A和哈密顿量之间的泊松括号运算也常用刘维尔算符表示 $$ \frac{dA}{dt}= \{A,H\}\equiv iLA $$ 该方程的形式解为 $$ A(t)=e^{iL}A(t=0) $$

可以证明,刘维尔算符是厄米算符,故上述指数算符是幺正算符。我们以后会大量使用刘维尔算符进行推导。

刘维尔定理

就像可以定义质量密度和电荷密度一样,也可以定义相空间的相点密度。首先定义相空间的体积元 \begin{equation} d \Gamma = dq_1 dq_2 \cdots dq_s dp_1 dp_2 \cdots dp_s = dq dp \end{equation} 体积元并不是数学上严格的微分,其大小的选取满足“宏观小”和“微观大”的要求。上式中第二个等号右边是常用的简写形式。如果在该体积元中相点的个数为 [math]dN[/math],那么就可以定义该体积元所在之处的相点密度: \begin{equation} \rho(q, p) = \frac{dN}{d \Gamma} \end{equation} 上式中的相点 [math](q, p)[/math] 可以是体积元 [math]d \Gamma[/math] 中的任意一点。如果将相空间类比于普通三维空间,将相点类比为三维空间中流体的质点,那么相点密度就对应于流体的质量密度。这个类比对下面的讨论是非常有用的。

下面我们来证明一个在统计物理中非常重要的定理:刘维尔定理: \begin{equation} \frac{d \rho}{dt} = 0. \end{equation} 这个定理是说:相点密度不随时间的变化而变化。首先,我们要清楚的是,虽然我们上面将相点密度写成了广义坐标和广义动量的函数,并不能说相点密度对时间的导数显然是零,因为广义坐标和广义动量可以是时间的函数。确切地讲,刘维尔定理说的是对任意的相轨迹, 相点密度是一个常量。

正如上面提到的,我们将相空间和其中的相点比作流体。对于这样的“流体”,我们有如下守恒定律: \begin{equation} \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial q_{\alpha}} (\rho \dot{q}_{\alpha}) + \frac{\partial}{\partial p_{\alpha}} (\rho \dot{p}_{\alpha}) = 0. \end{equation} 类似的公式在电磁学中叫做电荷守恒定律。因为相点是不会凭空产生和消失的,那么这个“流体”还不是一般的流体,而是不可压缩流体。对于不可压缩流体,上式中的“散度”处处为零: \begin{equation} \frac{\partial}{\partial q_{\alpha}} (\rho \dot{q}_{\alpha}) + \frac{\partial}{\partial p_{\alpha}} (\rho \dot{p}_{\alpha}) = 0. \end{equation} 于是,我们有 \begin{equation} \frac{\partial \rho}{\partial t} = 0. \end{equation} 这就是说,相点密度不显含时间。这也是我们将相点密度写成 [math]\rho(q, p)[/math],而不是 [math]\rho(q, p, t)[/math] 的原因。

要证明刘维尔定理,就是要证明下式: \begin{equation} \frac{d \rho}{dt} = \frac{\partial \rho}{\partial q_{\alpha}} \dot{q}_{\alpha} + \frac{\partial \rho}{\partial p_{\alpha}} \dot{p}_{\alpha} = 0. \end{equation} 将此式与上面“散度”为零的式子比较可知,只要能证明下式就大功告成了: \begin{equation} \frac{\partial \dot{q}_{\alpha} }{\partial q_{\alpha}} + \frac{\partial \dot{p}_{\alpha} }{\partial p_{\alpha}} = 0. \end{equation} 根据哈密顿正则方程,该式显然是恒成立的。刘维尔定理证毕。

习题

简谐振子的拉格朗日方程

  • 写出简谐振子的拉格朗日量,并由拉格朗日方程推导出其运动方程。

答案: \begin{equation} L = \frac{1}{2} m \dot{x}^2 - \frac{1}{2} k \dot{x}^2 \end{equation}


简谐振子的哈密顿方程

  • 题目:

从简谐振子的拉格朗日量推导出它的哈密顿量。该哈密顿量与简谐振子的机械能有什么关系?用哈密顿正则方程推导简谐振子的运动方程。

简谐振子的相空间

  • 题目:

简谐振子的相空间是什么样的,其中的一条相轨迹又是什么样的?

两条不同的相轨迹不可能有交点

  • 题目:

证明两条不同的相轨迹绝对不会相交于某一个相点。

热力学 (已修订)

基本概念

热力学系统

  • 一个热力学系统 (thermodynamic system) 是大量粒子的集合,常简称为系统 (system)。系统之外的所有物质称为环境 (environment)。
  • 如果系统与环境既不能交换能量也不能交换物质,则称其为一个孤立系统 (isolated system);若可交换能量但不能交换物质,则称为闭合系统 (closed system);若即可交换能量也可交换物质,则称为开放系统 (open system)。交换的能量可以是热能 (thermal energy),称为热量 (heat),也可以是机械能 (mechanical energy),称为 (work)。
  • 如果我们说一个系统达到了热力学平衡 (thermodynamic equilibrium),那么它的各个部分 (叫做子系统) 之间一定同时达到了
    • 热平衡 (thermal equilibrium)
    • 力学平衡 (mechanical equilibrium)
    • 扩散平衡 (diffusive equilibrium)

热力学第零定律和温度

  • 热力学中一共有四条基本经验定律,其中在时间上最后提出但在逻辑上排在第一个的叫热力学第零定律 (the zeroth law of thermodynamics)。 该定律说如果两个系统分别与第三个系统能达到热平衡,那么这两个系统也一定能达到热平衡。这说明处于热平衡的两个系统之间有一个相等的量。这个量就是热力学中最重要的物理量之一:温度 (temperature)。这就是温度的定义之一。
  • 通过规定两个特定系统的温度值,可以建立温标 (temperature scale)。一般选择一个大气压下水的冰点和沸点作为这两个特定系统。以下是两个常用温标:
    • 摄氏 (Celsius) 温标规定冰点的温度为 0 摄氏度,沸点的温度为 100 摄氏度,而且将这中间的温度均分为 100 份。这是我们生活中常用的温标。
    • 开尔文 (Kelvin) 温标规定冰点的温度为 273.15 K,沸点的温度为 373.15 K,而且将这中间的温度均分为 100 份。这是本书后面使用的温标。

状态与过程

  • 考虑一个简单的孤立或者闭合系统,达到热力学平衡时,可由三个物理量描述。它们是:压强 [math]p[/math]、体积 [math]V[/math] 和温度 [math]T[/math]。给定这三个量,就确定了系统的一个状态 (state)。事实上,这三个量并非完全独立,而是由状态方程 (equation of state) 联系着的。该方程可写为

$$ f(p, V, T) = 0. $$ 其中, [math]f[/math] 是一个特定的三元函数。当然我们可以将其中任何一个量写成另外两个量的两元函数,比如 [math]V=V(p,T)[/math]。这样的两元函数在由 [math]p[/math][math]V[/math][math]T[/math] 这三个参数构成的三维空间中可以表示为一个曲面。系统的某一个状态就对应于该曲面上的一个点。

  • 如果系统的状态发生了变化,我们称系统经历了一个过程 (process)。如果系统在其状态发生变化时始终无限接近平衡态,那么系统经历的过程为准静态过程 (quasi-static process)。显然,准静态过程对应于状态方程曲面上的一条曲线。

理想气体

  • 理想气体 (ideal gas) 是严格满足如下定律的系统:
    • 玻意耳-马略特 (Boyle-Mariotte) 定律:等温时压强反比于体积。
    • 查理 (Charles) 定律:等容时压强正比于温度。
    • 盖吕萨克 (Gay-Lussac) 定律:等压时体积正比于温度。
  • 通过这些定律以及标准状态 (压强为一个大气压,温度为 0 摄氏度) 下的摩尔体积 (约 [math]22.4\times10^{-3}[/math] m[math]^3[/math]/mol) 可以定出理想气体的状态方程

$$ pV = \nu R T. $$ 这里,[math]\nu[/math]物质的量 (amount of substance),[math]R \approx 8.31[/math] J/mol/K 是理想气体普适常量 (universal constant of ideal gas)。它与玻尔兹曼 (Boltzmann) 常量 [math]k_B \approx 1.38\times 10^{-23}[/math] J/K 由下式联系: \begin{equation} R \equiv N_A k_B. \end{equation} 这里,[math]N_A \approx 6.02\times 10^{23}[/math] mol[math]^{-1}[/math]阿伏伽德罗 (Avogadro) 常量。注意到 [math]\nu N_A =N[/math] 就是系统的粒子数,我们可将理想气体的状态方程改写成 \begin{equation} p V = N k_B T. \end{equation} 再定义 \begin{equation} n=N/V \end{equation} 为系统的数密度 (number density),我们又可以将理想气体的状态方程改写成 \begin{equation} p = n k_B T. \end{equation} 对于单原子理想气体,可以证明: \begin{equation} p = \frac{1}{3}n m \langle\vec{v}^2\rangle. \end{equation} 其中,[math]m[/math] 是一个原子的质量,[math]\langle\vec{v}^2\rangle[/math] 是原子速度平方的平均值。于是,我们有 \begin{equation} \frac{3}{2}k_BT = \frac{1}{2} m \langle\vec{v}^2\rangle. \end{equation} 这就是温度的微观意义之一:温度是原子平均平动动能的量度。

热力学第一定律

  • 热力学第一定律 (the first law of thermodynamics) 是说,在一个过程中,如果环境对系统做了功,也传了热,那么系统的内能 (internal energy) 的增加量 [math]\Delta E[/math] 就等于环境对系统做的功 [math]W[/math] 和传给系统的热 [math]Q[/math] 的和:

\begin{equation} \Delta E = Q + W. \end{equation} 当然,如果系统对环境做功,则约定 [math]W\lt 0[/math];如果系统传给环境热量,则约定 [math]Q\lt 0[/math]。这个定律其实就是能量守恒定律在热力学系统中的应用。

  • 如果考虑的是无限小过程,那么可以将热力学第一定律写成

\begin{equation} d E = \delta Q + \delta W. \end{equation} 功和热量都不是状态量,而是过程量,即它们依赖于具体的过程。或者说,功和热量不能描述系统的状态。这点与内能不同。从数学的角度来说,状态量的微分是恰当微分 (exact differential),而过程量的微分则不是。所以,为了区分,我们用 [math]dE[/math] 表示内能的微分,而用 [math]\delta W[/math][math]\delta Q[/math] 表示微小的功和热量。

热容

  • 系统在吸热时一般温度会升高 (有反例)。因为热量与过程有关,所以将一个系统的温度升高一定的值所需要的热量依赖于系统所经历的过程。指定一个过程,就可以定义一个叫做热容 (heat capacity) 的物理量:

\begin{equation} C = \frac{\delta Q}{d T}. \end{equation} 常见的两个过程是等容过程等压过程,对应的热容分别为等容热容等压热容

  • 如果体积固定,系统与外界互不做功,热力学第一定律告诉我们 [math]\delta Q = d E[/math]。从而,等容热容可表达为

\begin{equation} C_V =\left(\frac{~\delta Q}{dT}\right)_V = \left(\frac{\partial E}{\partial T}\right)_V. \end{equation}

  • 如果压强固定 (但体积不固定),系统要对环境做功 [math]pdV[/math],此时热力学第一定律告诉我们 [math]\delta Q = d E + p d V = d (E + p V)[/math]。从而,等压热容可表达为

\begin{equation} C_p = \left(\frac{~\delta Q}{dT}\right)_p = \left(\frac{\partial H}{\partial T}\right)_p. \end{equation} 其中,我定义了一个新的类似于内能的热力学函数, (enthalpy): \begin{equation} H = E + p V. \end{equation}

  • 综上可知:
    • 等容过程中系统吸收的热量等于其内能的增加量;
    • 等压过程中系统吸收的热量等于其焓的增加量。

绝热过程

  • 若是一个过程中,系统与环境没有热交换,那么该过程叫做绝热过程 (adiabaric process)。对理想气体来说,绝热过程可以由下式描述:

\begin{equation} p V^{\gamma} = \text{constant}. \end{equation} 其中, \begin{equation} \gamma \equiv \frac{C_p}{C_V}. \end{equation} 是等压热容和等容热容的比值,叫做绝热指数。

  • 以上结果的证明留作习题。

循环过程与卡诺 (Carnot) 循环

  • 循环过程 (cyclic process) 是系统终态等于初态的过程。在 [math]p-V[/math] 图中,循环过程对应于一个闭合路径。若闭合路径为顺时针方向,则系统对环境做净功并从环境中吸净热,对应于现实中的热机 (heat engine);反之,环境对系统做净功并从系统中吸净热,对应于现实中的制冷机 (refrigerator)。
  • 理论上最重要的循环过程为理想气体的卡诺 (Carnot) 循环 (这里需要一个图,但是我先不画)。图中的 [math]1\rightarrow 2\rightarrow 3\rightarrow 4\rightarrow 1[/math] 为正向 (顺时针) 卡诺循环,对应卡诺热机,而 [math]1\rightarrow 4\rightarrow 3\rightarrow 2\rightarrow 1[/math] 为逆向 (逆时针) 卡诺循环,对应卡诺制冷机。其中 [math]1 \rightarrow 2[/math][math]3 \rightarrow 4[/math] 都是等温过程,温度分别为较高温度 [math]T_{1}[/math] 和较低温度 [math] T_{2}[/math]。这两个温度由环境中的热浴 (heat bath) 来保持。另外两个部分 ([math] 2 \rightarrow 3[/math][math]4 \rightarrow 1[/math]) 都是绝热过程。
  • 根据热力学第一定律,对于卡诺热机,系统将在 [math]1 \rightarrow 2[/math] 的等温过程中从温度为 [math]T_{1}[/math] 的高温热源吸热,并在 [math]3 \rightarrow 4[/math] 的等温过程中向温度为 [math]T_{2}[/math] 的低温热源放热。计系统吸入和放出的热量为 [math]Q_{1}[/math][math]|Q_{2}|[/math] ([math]Q_{2} \lt 0[/math]),并定义热机的效率 (efficiency) [math]\eta[/math] 为系统所做净功与从高温热源所吸热量的比值:

\begin{equation} \label{equation:eta_1} \eta = \frac{ Q_{1} - |Q_{2}| } { Q_{1} }. \end{equation}

  • 很容易证明 (留作习题),这个效率其实只与温度有关而且总小于 1:

\begin{equation} \label{equation:eta_2} \eta = \frac{ T_{1} - T_{2} } { T_{1} } = 1- \frac{ T_{2} } { T_{1} }. \end{equation}

  • 为什么热机的效率总小于 1 呢?即:为什么系统不能把吸收的热量全部转化为功呢?这是热力学第一定律回答不了的问题。要回答这个问题,我们需要学习热力学第二定律。

热力学第二定律

可逆与不可逆过程

  • 有些过程不违背热力学第一定律 (能量守恒定律),但却不可能发生。比如,一个容器中的气体分子不会自发地全部聚集到容器的某一部分,即使这样做也不违背热力学第一定律。反过来,如果容器中的气体分子一开始是被限制在容器的某一部分,在限制解除后,它们却能够自发地充满整个容器。我们把气体分子自发地充满容器的过程称为自由膨胀过程。由上面的描述可知,这个过程不能自发地反过来进行。我们称这样的过程是不可逆的。在 19 世纪中叶,科学家们从这样的现象中总结了一个定律,叫热力学第二定律。根据这个定律,就能理解或者证明一个过程是否可逆。人们总结出的实际经验是:可逆过程必须是准静态 (即不要太快) 的而且是无耗散 (即不要有摩擦之类) 的;不同时满足这两个条件的过程就是不可逆过程

热力学第二定律的两种表述

  • 热力学第二定律可以有多种表示方法,但如下两种特别有名:
    • 开尔文表述: 不可能将从单一热源吸收的热量全部转变为功而不产生其它影响。即:第二类永动机不可能被制造出来。
    • 克劳修斯 (Clausius) 表述: 不能将热量从一个低温物体传给一个高温物体而不产生其它影响。即:热量总是自发地从高温物体传给低温物体。
  • 可以用反证法证明:这两种表述是等价的 (留作习题)。

卡诺定理

也可以用反证法证明 (留作习题):卡诺热机是所有热机中效率最高的。这叫做卡诺定理

克劳修斯 (Clausius) 定理

比较上面关于卡诺热机效率的两个公式可知 (注意 [math]Q_2\lt 0[/math]) \begin{equation} \frac{Q_1}{T_1} + \frac{Q_2}{T_2} = 0. \end{equation} 该式用微积分的语言可以写成 \begin{equation} \oint \frac{~d Q}{T} = 0. \end{equation} 克劳修斯证明了该式不光对卡诺循环成立,而且对任何可逆循环过程都成立。对不可逆循环,克劳修斯证明下式成立: \begin{equation} \oint \frac{~d Q}{T} < 0. \end{equation} 以上两式分别称为克劳修斯等式和不等式。这两个式子表述的内容被称为克劳修斯定理

  • 克劳修斯定理的证明留作习题。

根据微积分知识,适用于可逆过程的克劳修斯等式意味着 [math]\frac{\delta Q}{T}[/math] 是一个恰当微分。或者说,积分 [math]\int_A^B \frac{\delta Q}{T}[/math] 与从初态 [math]A[/math] 到终态 [math]B[/math] 的路径无关。一个恰当微分总可以写成某个函数 (或者说变量;变量与函数的概念是相对的) 的全微分。我们用 [math]dS[/math] 来表示这个全微分: \begin{equation} dS \equiv \frac{\delta Q}{T}. \end{equation} 中国老辈物理学家为这个新的热力学函数翻译了个不错的名字: (entropy)。我们可以这样理解这个名字:第一,它是一个“商”,即热量除以温度;第二,它有个“火”字旁,表明它与热量有关。

  • 用这个记号,就可以把热力学第一定律写成一个只涉及热力学变量的形式:

\begin{equation} dE = TdS - pdV. \end{equation} 这个式子在热力学中特别重要,常被称为热力学基本方程。对热力学理论的探索常常从这个方程出发。

  • 从这个方程可以看出,熵是一个广延量,或则说可加量。更重要的是,由于熵是一个状态量,除了一个积分常数 (类似于势能的零点) 它的值不依赖于系统所经历的过程。
  • 克劳修斯定理 (即克劳修斯等式和不等式) 可以说是热力学第二定律的一种数学表述。然而,该数学表述是基于一个循环过程的。很多情况下,考虑一个普通的过程或者一个无限小过程是更方便的。通过适用于可逆过程的克劳修斯等式,我们确立了一个恰当微分 [math]dS=\frac{\delta Q}{T}[/math]。对此积分,便可得对应于一个从态 [math]A[/math] 到态 [math]B[/math] 的有限的可逆过程中的熵变:

\begin{equation} S_B-S_A = \int_{A}^{B} \frac{\delta Q}{T}. \end{equation} 现在呢,我们假设通过一个任意的 (可逆的或者不可逆的) 过程让系统从态 [math]B[/math] 回到态 [math]A[/math]。根据克劳修斯定理,我们有 \begin{equation} S_B-S_A + \int_{B}^{A} \frac{\delta Q}{T} = \int_{A}^{B} \frac{\delta Q}{T} + \int_{B}^{A} \frac{\delta Q}{T} = \oint \frac{\delta Q}{T} \leq 0, \end{equation} 即 \begin{equation} S_A-S_B \geq \int_{B}^{A} \frac{\delta Q}{T}. \end{equation} 令 [math]A[/math] 状态无限接近 [math]B[/math] 状态,即得 \begin{equation} dS \geq \frac{\delta Q}{T}. \end{equation} 从推导的过程可知,这里的等号和不等号分别对应于克劳修斯定理中的等号和不等号,从而分别对应于可逆与不可逆过程。由于与克劳修斯定理等价,这个不等式也可以被当做热力学第二定律的数学表述

  • 特别地,如果过程是绝热的,即 [math]\delta Q=0[/math],那么有

\begin{equation} dS \geq 0. \end{equation} 该式常被称为熵增加原理。一个更特殊的情况是孤立系统中的过程。由于孤立系统不与环境交换物质与能量 (当然,包括热能),无论它的内部发生何种变化,它的熵只能增加或不变,不能减少。

在结束这一节之前,我们来重新考察一下在本节开头提到过的一个不可逆过程:气体的自由膨胀。为简单起见,我们考虑单原子理想气体的绝热自由膨胀 (即假定理想气体被装在绝热容器里面)。如果气体的体积从 [math]V_1[/math] 自由地膨胀至 [math]V_2[/math],那么,由于其温度不变 (气体与环境无热交换,自由膨胀的理想气体与环境也无功的交换,故内能不变,从而温度不变),系统熵的增量为 ([math]N[/math] 为气体分子数) \begin{equation} \Delta S = Nk_B \ln \frac{V_2}{V_1} > 0. \end{equation} 于是,该绝热过程是不可逆的,与我们的直觉一致。这里,我们用到了理想气体熵增加量的公式 (见习题)。

热力学函数和关系

开放系统和化学势

  • 到目前为止,我们只研究了与环境没有物质交换的孤立和闭合系统。对于开放系统,其粒子数 [math]N[/math] 是允许变化的。为了能处理粒子数变化的过程,我们将 [math]N[/math] 也视作一个热力学变量并将热力学基本方程推广为

\begin{equation} \boxed{dE = T dS - P dV + \mu dN}. \end{equation}

这里,我们引入了一个新的强度量[math]\mu[/math],叫做化学势 (chemical potential)。从上述内能的全微分可以看出, \begin{equation} \mu = \left( \frac{\partial E}{\partial N} \right)_{S, V}. \end{equation} 这个式子可以看作是化学势的一个定义。可以证明,理想气体的化学势是负的。证明留作练习。

欧拉 (Euler) 方程

热力学基本方程告诉我们内能是熵、体积以及粒子数的一个函数: \begin{equation} E=E(S, V, N). \end{equation} 该式中所有的变量都是广延量,意味着内能是这些变量的齐次函数 (homogeneous function),即 \begin{equation} E(\lambda S, \lambda V, \lambda N) = \lambda E(S, V, N). \end{equation} 其中,[math]\lambda[/math] 是一个任意的参数。将该式对参数 [math]\lambda[/math] 求导可得 \begin{align} &\frac {\partial E(\lambda S, \lambda V, \lambda N)}{ \partial(\lambda S)} \frac {\partial (\lambda S)} {\partial \lambda} + \frac {\partial E(\lambda S, \lambda V, \lambda N)}{\partial (\lambda V)} \frac {\partial (\lambda V)} {\partial \lambda} + \frac {\partial E(\lambda S, \lambda V, \lambda N)}{\partial (\lambda N)} \frac {\partial (\lambda N)} {\partial \lambda} \nonumber \\ &= E(S, V, N). \end{align} 再将 [math]\lambda[/math] 取为 1,可得 \begin{equation} E = \frac {\partial E}{\partial S} S + \frac {\partial E}{ \partial V} V + \frac {\partial E}{\partial N} N. \end{equation} 另外,从推广的热力学基本方程可知如下关系: $$ T = \left( \frac{\partial E}{\partial S} \right)_{V, N} $$ $$ p = -\left( \frac{\partial E}{\partial V} \right)_{S, N} $$ 以及 $$ \mu = \left( \frac{\partial E}{\partial N} \right)_{S, V} $$ 于是,我们得到一个重要的方程: \begin{equation} \boxed{E = TS - pV + \mu N}. \end{equation} 该方程被称为欧拉 (Euler) 方程

吉布斯-杜安 (Gibbs-Duhem) 关系

由欧拉方程可以推导另一个重要的关系式。对欧拉方程两边作全微分,可得 \begin{equation} dE = T dS + SdT - p dV -V dp + \mu dN + N d \mu. \end{equation} 将此与推广的热力学基本方程对照便知 \begin{equation} \boxed{SdT - V dp + N d \mu = 0}. \end{equation} 该方程被称为吉布斯-杜安关系。这个关系告诉我们:三个强度量 (压强,温度,化学势) 中只有两个是独立的。因此,我们可以说:这样的系统的热力学自由度为 2。

从能量函数到熵函数

根据欧拉方程,我们可以将熵表达为能量、体积以及粒子数的函数: \begin{equation} S = S(E, V, N) = \frac{E+pV-\mu N}{T}. \end{equation} 其全微分可由热力学基本方程得到: \begin{equation} dS = \frac{1}{T} dE + \frac{p}{T} dV - \frac{\mu}{T} dN. \end{equation} 于是,我们从能量函数过渡到了熵函数。从熵函数的全微分出发,可以更方便地讨论热力学第二定律,因为热力学第二定律更多地是与熵,而不是能量相联系。

平衡条件

  • 我们在本讲第一节就提到了三个热力学平衡:热平衡,力平衡,以及扩散平衡。热力学第零定律告诉我们一个系统达到热平衡时内部温度处处相等,而基于力学的直觉告诉我们系统达到力平衡时应该处处压强相等。在这一小节,我们将从热力学第二定律的数学表述之一,即熵增加原理来更为严格地推导热平衡与力平衡的条件,以及我们还没有讨论过的扩散平衡条件。
  • 一个孤立系统是不受外界干扰的,无论它的初始条件如何,在等待足够长的时间之后,它总会达到一个固定的状态,宏观上表现为达到热力学平衡。这与上一节得到的孤立系统的熵总是趋于最大化的结论不谋而合。所以,对热力学平衡的向往是孤立系统中熵增加的动力。如果孤立系统的熵在增加,那么它就还未达到热力学平衡;如果孤立系统的熵达到了一个最大值,不再增加了,它就达到了热力学平衡;在达到热力学平衡之后,孤立系统的熵就会保持最大值,不再变化。
  • 为简单起见,我们考虑一个孤立系统,并用一个假想的边界将该系统分为两个子系统:[math]A[/math][math]B[/math]。显然,这个两个子系统都是开放系统,它们之间能够交换物质与能量。我们记 [math]A[/math] 系统中的热力学量为 [math]S_A[/math][math]E_A[/math][math]V_A[/math][math]N_A[/math][math]T_A[/math][math]p_A[/math]、和 [math]\mu_A[/math];记 [math]B[/math] 系统中的热力学量为[math]S_B[/math][math]E_B[/math][math]V_B[/math][math]N_B[/math][math]T_B[/math][math]p_B[/math]、和 [math]\mu_B[/math]。 由于熵是可加量 (广延量),整个孤立系统的熵的全微分可写为

\begin{equation} d S = d S_A+d S_B = \frac{1}{T_A} dE_A + \frac{p_A}{T_A} dV_A - \frac{\mu_A}{T_A} dN_A + \frac{1}{T_B} dE_B + \frac{p_B}{T_B} dV_B - \frac{\mu_B}{T_B} dN_B. \end{equation} 因为整个系统的粒子数、体积、以及内能都是守恒的,即 [math]dN_A+dN_B=0[/math][math]dV_A+dV_B=0[/math] 以及 [math]dE_A+dE_B=0[/math],我们有 \begin{equation} d S = \left(\frac{1}{T_A} - \frac{1}{T_B} \right) dE_A + \left(\frac{p_A}{T_A} - \frac{p_B}{T_B} \right) dV_A - \left(\frac{\mu_A}{T_A} - \frac{\mu_B}{T_B} \right) dN_A. \end{equation}

  • 若整个孤立系统达到了热力学平衡,则有 [math]dS=0[/math],从而有如下平衡条件:

\begin{equation} T_A = T_B \quad (\text{热平衡}), \end{equation} \begin{equation} p_A = p_B \quad (\text{力平衡}), \end{equation} \begin{equation} \mu_A = \mu_B \quad (\text{扩散平衡}). \end{equation} 前两个平衡条件与我们的之前的结果是一致的。最后一个平衡条件是说孤立系统在达到扩散平衡时,其内部的化学势处处相等。

  • 如果系统仍处在向热力学平衡奔波的路上,熵增加原理,即 [math]dS\gt 0[/math] 告诉我们整个系统的广延量 (即内能、体积、以及粒子数) 会在子系统之间按照如下规则重新分配:

\begin{equation} T_A > T_B \Longrightarrow d E_A < 0, \end{equation} \begin{equation} T_A = T_B ~\text{且}~ p_A > p_B \Longrightarrow d V_A > 0, \end{equation} \begin{equation} T_A = T_B ~\text{且}~ \mu_A > \mu_B \Longrightarrow d N_A < 0. \end{equation} 也就是说,在趋向平衡的过程中,具有较高温度的子系统会失去内能以降低温度,具有较高压强的子系统会扩展体积以降低压强,具有较高化学势的子系统会失去物质 (粒子) 以降低化学势。所以,化学势是对粒子扩散趋势的一种量度

热力学势

  • 无论是能量的全微分还是熵的全微分,其独立变量 (也叫自然变量) 都是广延量。然而,实验上更容易控制的是强度量。有没有一种方法将部分或全部独立变量用强度量替代呢?勒让德变换可以帮助我们做到这一点。勒让德变换可以将一个函数变换为另外相关的一个函数。具体地,新的函数是将原来的函数减去该函数的一个独立变量与对应的偏导数的乘积而得到的。
  • 我们的出发点是欧拉公式:

\begin{equation} E = E(S, V, N) = TS -pV + \mu N \end{equation} 和它的全微分,即热力学基本方程: \begin{equation} dE = T dS - p dV + \mu dN. \end{equation}

  • 考虑能量函数 [math]E = E(S, V, N)[/math] 并取 [math]V[/math] 为需要变换的独立变量,对应的勒让德变换将能量函数变换为焓函数:

\begin{equation} H = H(S,p,N)=E - \left(\frac{\partial E}{\partial V}\right)_{S,N} V= E + p V = TS + \mu N. \end{equation} 其中,我们在最后一步用了欧拉公式。你可能也注意到上式暗示了焓函数的独立变量是 [math]S[/math][math]p[/math][math]N[/math]。这一论断能由下式证实 (第二个等号运用了热力学基本方程): \begin{equation} dH = dE + pdV + Vdp = T dS + V dp + \mu dN. \end{equation} 对焓的表达式 [math]H=TS+\mu N[/math] 求微分并运用吉布斯-杜安关系也能得到这个结果。

亥姆霍兹 (Helmholtz) 自由能

  • [math]S[/math] 为需要变换的独立变量,对能量函数进行勒让德变换,就可以得到亥姆霍兹 (Helmholtz) 自由能函数:

\begin{equation} F = F(T,V,N)=E - \left(\frac{\partial E}{\partial S}\right)_{V,N} S= E - T S = -pV + \mu N. \end{equation} 它的独立变量为 [math]T[/math][math]V[/math][math]N[/math],因为 \begin{equation} dF = dE - TdS - S dT = -S dT - pdV + \mu dN. \end{equation}

吉布斯 (Gibbs) 函数

如果同时针对 [math]V[/math][math]S[/math] 两个变量进行勒让德变换,则可得到吉布斯 (Gibbs) 函数: \begin{equation} G = G(T, p, N) = E - T S + p V = \mu N. \end{equation} 它的独立变量为 [math]T[/math][math]p[/math][math]N[/math],因为 \begin{equation} dG = - S d T + V dp + \mu dN. \end{equation} 其实,吉布斯函数也可以认为是从焓对熵进行一次勒让德变换得到的:[math]G=H-TS[/math],或者从亥姆霍兹自由能对体积进行一次勒让德变换得到的:[math]G=F+pV[/math]。不管怎样看待,吉布斯函数的全微分都有上式中的形式。

巨正则函数

最后,如果取 [math]N[/math] 为需要变换的独立变量并对亥姆霍兹自由能进行勒让德变换,我们将得到 巨正则函数 (Grand canonical function): \begin{equation} \Phi = \Phi(T, V, \mu) = F - \mu N = - p V. \end{equation} 它的独立变量为 [math]T[/math][math]V[/math][math]\mu[/math],因为 \begin{equation} d\Phi = - S d T - p d V - N d \mu. \end{equation}

麦克斯韦 (Maxwell) 关系

  • 内能和上述四个通过勒让德变换得到的函数统称为热力学势 (thermodynamic potentials)。 我们能用勒让德变换定义更多的函数,但我们以后只会用到这几个。
  • 因为热力学势都是系统的状态变量,它们的全微分都是恰当的。这意味着它们对任何两个独立变量的二阶偏导数与求偏导的次序无关。例如,根据

$$dE = T dS - p dV + \mu dN$$ 我们有 \begin{equation} \frac{\partial}{\partial V} \left(\frac{\partial E}{\partial S}\right) = \frac{\partial}{\partial S} \left(\frac{\partial E}{\partial V}\right). \end{equation} 将等式两边括号内的一阶偏导数用合适的变量替代,可得 \begin{equation} \left( \frac{\partial T}{\partial V} \right)_{S, N} = - \left( \frac{\partial p}{\partial S} \right)_{V, N}. \end{equation} 这样得到的关系称为麦克斯韦 (Maxwell) 关系。显然,可以有很多麦克斯韦关系,完全的推导留作习题。

最小能量原理

  • 前面,我们根据熵增加原理讨论了孤立系统趋于平衡的过程。然而,孤立系统的模型往往不是讨论一个特定问题时的最佳选择。这里,我们来考察其它类型的系统 (闭合系统和开放系统) 趋于平衡的过程。
  • 一个常见的例子是体积和粒子数固定的闭合系统,它与一个热浴接触而保持一个恒定的温度。 对这样的非孤立系统,熵增加原理不再适用,因为系统与环境 (主要是热源) 可能会交换热量。为了研究此时系统的热力学演化行为,我们必须从热力学第二定律的更为一般的表达式,即 [math]TdS \geq \delta Q[/math] 出发。一方面,由于温度是恒定的,我们有 [math] TdS = d(TS)[/math];另一方面,由于粒子数和体积固定,热力学第一定律告诉我们 [math]\delta Q = dE[/math]。结合这三个式子,并注意到亥姆霍兹自由能的定义 ([math]F=E-TS[/math]),我们得到如下不等式:

\begin{equation} d F \leq 0. \end{equation} 该式被称为自由能最小原理。它是说,一个粒子数、体积和温度固定的系统总是朝着自由能减小的方向演化的。如果系统的自由能还未达到最小值,那么它还未达到热力学平衡;如果系统的自由能达到了最小值,那么它就达到了热力学平衡。

  • 另一个常见的例子是粒子数、压强和温度固定的系统。类似地,我们可以根据热力学第一和第二定律得到不等式 [math]d (TS) \geq dE + d(pV)[/math]。利用吉布斯函数的定义可以将该不等式表达为

\begin{equation} d G \leq 0. \end{equation} 这就是吉布斯函数最小原理,即一个粒子数、压强和温度固定的系统总是朝着吉布斯函数减小的方向演化。吉布斯函数减小时体系还未达到热力学平衡,吉布斯函数不变时体系便达到了热力学平衡。

  • 以上两个原理统称为最小能量原理,它是力学中的最小势能原理在热力学中的推广。由于我们是从热力学第一和第二定律出发推导该原理的,它也可以被当做是热力学第二定律的数学表述之一。它与熵增加原理在本质上是等价的,而且在应用中是互为补充的。

习题

理想气体绝热过程方程

  • 题目:请推导理想气体绝热过程方程:

\begin{equation} p V^{\gamma} = \text{constant}. \end{equation} 其中, \begin{equation} \gamma \equiv \frac{C_p}{C_V}. \end{equation} 是等压热容和等容热容的比值,叫做绝热指数。

卡诺热机效率

题目:请推导正文中给出的理想气体卡诺热机的效率公式。

理想气体的熵

  • 题目:

考虑单原子理想气体,其内能为体系的动能 $$E=\frac{3}{2}Nk_BT.$$ 证明从温度为 [math]T_1[/math]、体积为 [math]V_1[/math] 的状态变化到温度为 [math]T_2[/math]、体积为 [math]V_2[/math] 的状态的过程中,系统的熵增加量为 \begin{equation} S_2 - S_1 =Nk_B \ln \left[ \frac{V_2}{V_1} \left(\frac{T_2}{T_1}\right)^{3/2} \right]. \end{equation}

理想气体的化学式

  • 题目:证明理想气体的化学式为负。
  • 解答:根据上一题的结果,熵是随着 [math]N[/math][math]E[/math][math]V[/math] 这些变量的增大而增大的。因此,如果要在增大 [math]N[/math] 的同时将 [math]S[/math][math]V[/math] 固定,必须将 [math]E[/math] 减小。因此,根据化学势的定义式,理想气体的化学势必然是负的。

麦克斯韦关系式

请根据正文中讨论过的热力学势推导所有的麦克斯韦关系式。

经典统计力学 (已修订)

微正则系综理论

宏观态与微观态

  • 我们考虑一个孤立系统。这个系统具有固定的粒子数 [math]N[/math] (系统与外界无粒子交换)、体积 [math]V[/math](系统与外界无体积功的交换)以及能量 [math]E[/math](系统与外界也不交换热量,故能量不变)。为了简化讨论,我们还假设系统中的粒子之间没有相互作用,即考虑理想气体系统。还是为了讨论方便,我们假设粒子的能量可以取一系列离散的值 [math]E_i~(i=1, 2, \cdots)[/math],并设能量等于 [math]E_i[/math] 的粒子的数目为 [math]n_i[/math]。于是,我们有:

\begin{equation} N = \sum_i n_i, \end{equation} \begin{equation} E = \sum_i n_i E_i. \end{equation}

我们知道,从热力学的角度来看,系统的状态就由 [math]N[/math][math]V[/math][math]E[/math] 来描述。这些量叫做宏观量,而一组宏观量就确定了一个宏观态 (macroscopic state)。

  • 然而,从微观的角度来说,即使一个系统处于一个确定的宏观态,系统中各个粒子的运动状态也可能是不确定的。系统中所有粒子的运动状态的组合构成一个微观态 (microscopic state)。一个 [math]N[/math][math]V[/math][math]E[/math] 确定的宏观态可能具有多个微观态,而且微观态的个数一般来说是 [math]N[/math][math]V[/math][math]E[/math] 的函数。我们将这个函数记为

\begin{equation} \Omega = \Omega(N, V, E). \end{equation}

  • 上面只考虑了一个孤立系统。现在考虑一个由子系统 1 和子系统 2 构成的复合系统。令子系统 1 的宏观态由粒子数 [math]N_1[/math]、体积 [math]V_1[/math] 和能量 [math]E_1[/math] 表示,子系统 2 的宏观态由粒子数 [math]N_2[/math]、体积 [math]V_2[/math] 和能量 [math]E_2[/math] 表示。子系统 1 和 2 的微观状态数分别记为 [math]\Omega_1 = \Omega_1(N_1, V_1, E_1)[/math][math]\Omega_2 = \Omega_2(N_2, V_2, E_2)[/math]。设这两子个系统之间可以相互传热,但不能相互传递粒子,也不能相互做体积功。于是,每子个系统的粒子数和体积都不变,但能量可能随时间的推移而改变。设两个子系统之间的相互作用能可以忽略不计(为什么可以忽略不计?什么情况下不可以忽略不计?),则复合系统的总能量 [math]E[/math] 就是每个子系统的能量的和:

\begin{equation} E_1 + E_2 = E. \end{equation} 我们假定复合系统是孤立的,其总能量 [math]E[/math] 是一个常数。

  • 因为子系统 1 有 [math]\Omega_1(N_1, V_1, E_1)[/math] 个微观状态,子系统 2 有 [math]\Omega_2(N_2, V_2, E_2)[/math] 个微观状态,而子系统 1 的任何一个微观状态与子系统 2 的任何一个微观状态一起构成了复合系统的一个可能的微观状态,故由乘法原理可知,复合系统有

\begin{equation} \Omega(N_1, V_1, E_1, N_2, V_2, E_2) = \Omega_1(N_1, V_1, E_1)\Omega_2(N_2, V_2, E_2) \end{equation} 个微观状态。我们现在的问题是:总能量 [math]E[/math] 如何在两个子系统之间分配,才能让两个子系统之间达到热力学平衡?也就是说,如果一开始两个子系统之间没有达到热力学平衡,那么能量 [math]E_1[/math](或者能量 [math]E_2[/math])要如何改变才能使得两个子系统趋向热力学平衡?

等概率原理

  • 要回答上面的问题,必须对宏观态与微观态的对应作出一个假设。这个假设就是等概率原理。这个原理是说,一个孤立系统所有可能的微观态出现的概率都是相等的。等概率原理是统计力学中唯一的原理。这个原理正确与否只有实验能够判决。迄今为止的所有科学事实证明,这个原理是非常合理的。
  • 我们对等概率原理并不陌生。例如,当你说抛掷一个色子得到一点的概率是 [math]1/6[/math] 时你就运用了等概率原理。
  • 等概率原理如何运用到统计力学和热力学中呢?我们知道,一个系统在一定的宏观约束下具有多个可能的宏观态,每个宏观态对应一定的微观状态数。既然每一个微观状态——不管它属于哪一个宏观态——出现的概率都相等,那么那个具有最大微观状态数的宏观态出现的概率自然是最大了。这个具有最大微观状态数的宏观态叫做最可几态 (most probable state)。所以,随着时间的推移,系统将演化到最可几态。这个最可几态自然就是平衡态

熵的微观解释

  • 根据上面的等概率原理,当复合系统达到热力学平衡时,微观状态数 [math]\Omega(N_1, V_1, E_1, N_2, V_2, E_2)[/math] 应该具有最大值。因为这个问题中唯一可以变化的是 [math]E_1[/math] (或者 [math]E_2[/math]),这个最大值发生的条件是:

\begin{equation} \frac{\partial }{\partial E_1}\Omega(N_1, V_1, E_1, N_2, V_2, E_2) = \frac{\partial }{\partial E_1} [\Omega_1(N_1, V_1, E_1)\Omega_2(N_2, V_2, E_2)]=0. \end{equation} 在我们的问题中,[math]N_1[/math][math]N_2[/math][math]V_1[/math][math]V_2[/math] 都是常数,故可将上式简写为 \begin{equation} \frac{\partial }{\partial E_1} [\Omega_1(E_1)\Omega_2(E_2)]=0. \end{equation} 注意到 [math]E_1 + E_2 = E[/math],我们有 \begin{equation} \left( \frac{\partial \Omega_1(E_1)}{\partial E_1} \right) \Omega_2(E_2) - \Omega_1(E_1) \left( \frac{\partial \Omega_2(E_2)}{\partial E_2} \right) = 0, \end{equation} 即 \begin{equation} \frac{1}{\Omega_1(E_1)} \left( \frac{\partial \Omega_1(E_1)}{\partial E_1} \right)= \frac{1}{\Omega_2(E_2)} \left( \frac{\partial \Omega_2(E_2)}{\partial E_2} \right), \end{equation} 亦即 \begin{equation} \left( \frac{\partial \ln [\Omega_1(E_1)]} {\partial E_1} \right)_{N_1V_1}= \left( \frac{\partial \ln [\Omega_2(E_2)]} {\partial E_2} \right)_{N_2V_2}. \end{equation} 注意:我们在两边的求导操作中加上了之前省略不写的条件。

  • 上式就是两个子系统达到热力学平衡的条件。因为子系统之间不能传粒子也不能做功,这个条件具体地说就是热平衡的条件。这不就是热力学第零定律吗?因此,我们从统计力学推导出了热力学第零定律:两个子系统达到热平衡时必有某个量相等,而这个量就是温度。
  • 热平衡时子系统之间的温度相等,即 [math]T_1=T_2[/math],或者 [math]1/T_1=1/T_2[/math]。运用热力学基本方程,我有

\begin{equation} \left( \frac{\partial S_1} {\partial E_1} \right)_{N_1V_1} = \left( \frac{\partial S_2} {\partial E_2} \right)_{N_2V_2}. \end{equation} 比较以上两式,完全可以作出如下猜测: \begin{equation} S \propto \ln [\Omega(E)]. \end{equation} 也就是说,一个系统的某个宏观态的熵正比于该宏观态的微观状态数的对数。这里,我们去掉了表示系统的下标,因为这个公式对任何系统都是成立的。这就是熵的微观解释,是玻尔兹曼 (Boltzmann) 最先得到这个关系的。

如何确定上式中的比例常数呢?我们注意到玻尔兹曼常数 [math]k_B[/math] 与熵具有相同的量纲,而微观状态数的对数的量纲是 1,故可猜测这个比例常数就是 [math]k_B[/math]: \begin{equation} \boxed{S = k_B \ln [\Omega(E)]}. \end{equation} 历史上,是普朗克 (Planck) 最先写出这个等式的。这正是刻在玻尔兹曼墓碑上的公式。

  • 下面,我们以理想气体为例来验证这个比例常数确实就是玻尔兹曼常数。首先,可以确定,理想气体的微观状态数与气体的体积有如下关系:

\begin{equation} \Omega(N, V, E) \propto V^N. \end{equation} 再由热力学关系 ([math]p[/math] 是系统的压强;[math]T[/math] 是系统的温度) $$ \frac{p}{T} = \left( \frac{\partial S}{\partial V} \right)_{NE} $$ 可得 \begin{equation} pV = N k_B T. \end{equation} 这正是理想气体状态方程。于是可以确定,[math]k_B[/math] 就是以前在热力学中引入的玻尔兹曼常数。

  • 最后,我们注意到,熵的微观解释和等概率原理也包含了热力学第二和第三定律:
    • 根据等概率原理,一个孤立系统总是向微观状态数最大的宏观态演化。这就是用来表述热力学第二定律的熵增加原理。
    • 根据熵的公式 [math]S=k_B \ln \Omega[/math],熵有一个绝对的最小值 0,当且仅当 [math]\Omega=1[/math] 时取得。这为热力学第三定律提供了一个理论基础。

正则系综理论

正则系综的概率函数

上一讲研究的具有固定的粒子数、体积和能量的孤立系统对应于微正则系综中的一个系统。一个系综就是由若干个假想的与所研究的系统具有相同的宏观态的系统的集合。微正则系综在很多情况下都不方便使用。本讲研究一个更有用的系综:正则系综(canonical ensemble)。一听名字就知道正则系综更加正统,有用。

正则系综考虑的是 [math]M[/math] 个相同的具有一定粒子数、体积和温度的热力学系统。可以想象 [math]M[/math] 个相同的系统排成一个圈;相邻的系统之间有微弱的相互作用,使得所有的系统最终能处于同一温度。用 [math]M_i[/math] 表示在微观状态 [math]i[/math] 上的系统数目,[math]E_i[/math] 表示第 [math]i[/math] 个态的能量(需要假设能级不简并吗?我不是很清楚)。则总系统数和总能量为 \begin{equation} M = \sum M_i. \end{equation} \begin{equation} E_M = \sum M_i E_i. \end{equation} 知道了分布 [math]\{M_i\}[/math],并没有确定各个系统的状态。可以证明 (留作习题):对于一个给定的分布 [math]\{M_i\}[/math],系综的微观状态数为 \begin{equation} \Omega = \frac{M!}{\prod_i M_i!}. \end{equation} 从而,这个系综(作为一个很大的孤立系统)的熵为 \begin{equation} S_M = k_B \ln \Omega = k_B \ln \left(\frac{M!}{\prod_i M_i!}\right). \end{equation}

下面我们要问:哪个分布 [math]\{M_i\}[/math] 的概率最大?那个具有最大概率的分布就对应于平衡态。根据等概率原理,肯定是微观状态数最大的分布概率最大。因为微观状态数往往是巨大的,所以一般不直接求 [math]\Omega[/math] 的最大值,而是求 [math]\ln \Omega[/math] 的最大值。

对于一个多元函数 [math]\Omega(M_1, M_2, \cdots)[/math],要使其取最大值,自然的想法是令其对所有变量的导数都等于零: \begin{equation} \frac{\partial \ln \Omega}{\partial M_i} = 0. \end{equation} 然而,这个式子是错的,因为我们忽略了上面的约束条件 [math]M = \sum M_i[/math][math]E_M = \sum M_i E_i[/math]。有约束条件的极值问题一般用拉格朗日乘子法来解决。引入两个拉格朗日乘子 [math]\alpha[/math][math]\beta[/math]。极值条件可以写成 \begin{equation} \frac{\partial \ln \Omega}{\partial M_i} - \alpha \frac{\partial \sum_j M_j}{\partial M_i} - \beta \frac{\partial \sum_j M_j E_j}{\partial M_i} = 0. \end{equation}

在进一步推导之前,我们先证明一个Stirling公式,即当 [math]M \gg 1[/math] 时,有 \begin{equation} \ln M! = M(\ln M - 1) + O(\ln M). \end{equation} 其中,[math]O(\ln M)[/math] 代表与 [math]\ln M[/math] 同数量级的项。首先,我们计算积分 \begin{equation} I = \int_1^M \ln(x)dx = M(\ln M - 1) + 1. \end{equation} 可以证明: \begin{equation} \ln M! - \ln M < I < \ln M!. \end{equation} 于是有 \begin{equation} M(\ln M - 1) + 1 < \ln M! < M(\ln M - 1) + 1 + \ln M. \end{equation}

[math]M[/math] 趋近于无穷大时,所有的 [math]M_i[/math] 也趋近于无穷大。利用斯特林公式,有 \begin{equation} \ln \Omega \approx M(\ln M - 1)- \sum_i \left[ M_i(\ln M_i - 1) \right]. \end{equation}

可以证明 \begin{equation} \frac{\partial M}{\partial M_i} = 1. \end{equation} \begin{equation} \frac{\partial E_M}{\partial M_i} = E_i. \end{equation} \begin{equation} \frac{\partial \ln \Omega}{\partial M_i} = - \ln M_i. \end{equation} 于是,我们可以将极值条件写成 \begin{equation} \ln M_i = -\alpha - \beta E_i. \end{equation} 对上式两边取指数,可得 \begin{equation} M_i = \exp[-\alpha - \beta E_i]. \end{equation}

既然 [math]M_i[/math] 是处于微观态 [math]i[/math] 的系统的个数,那么自然地,一个系统处于微观态 [math]i[/math] 的概率为 \begin{equation} w_i = \frac{M_i}{M} = \frac{\exp[- \beta E_i]}{\sum_i \exp[-\beta E_i]}. \end{equation} 可见,常数 [math]\alpha[/math] 没有什么物理意义,在求概率的时候就被消掉了。上式中的分母叫做正则配分函数(canonical partition function),记为 [math]Z[/math]: \begin{equation} \boxed{Z = \sum_i \exp[-\beta E_i]}. \end{equation}

表面上看,配分函数只不过是一个归一化因子而已。然而实际上,配分函数包含了体系所有的热力学性质!后面马上会验证此论断。


正则系综中的热力学量

凭猜测引入温度

  • 那么,常数 [math]\beta[/math] 有什么物理意义呢?一方面,我们知道,正则系综中的温度是一个常数,所以可以猜测 [math]\beta[/math] 与温度有关。另一方面,从量纲的角度来看,可以进一步猜测:

\begin{equation} \boxed{\beta = \frac{1}{k_B T}}. \end{equation} 于是,正则配分函数可以写成 \begin{equation} \boxed{Z = \sum_i \exp \left[-\frac{ E_i } { k_B T} \right]}. \end{equation} 怎么确定这里引入的 [math]T[/math] 就是绝对温度呢?等到讲完后面的位力定理时我们就会明白。在那之前,先假设 [math]T[/math] 就是绝对温度,并考察几个热力学量的计算。

  • 可以证明,整个正则系综的熵可以化为如下形式:

\begin{equation} S_M = -k_B M \sum_i w_i \ln w_i. \end{equation} 于是,由熵的可加性得到单个系统的熵为 \begin{equation} S = -k_B \sum_i w_i \ln w_i. \end{equation} 这个公式叫做熵的吉布斯公式。

亥姆霍兹自由能

  • 将概率函数 [math]w_i[/math] 的表达式代入熵的吉布斯公式可得

\begin{equation} S = k_B \ln Z + \frac{E}{T}. \end{equation} 其中, \begin{equation} E = \sum_i w_i E_i \end{equation} 是系统能量的平均值。于是, \begin{equation} E - TS = -k_B T\ln Z. \end{equation} 上式右边就是亥姆霍兹自由能: \begin{equation} \boxed{F = -k_B T \ln Z}. \end{equation} 这样就将正则配分函数与系统的亥姆霍兹自由能联系起来了。这个联系的意义是深远的,因为我们知道,在粒子数、体积、温度固定的系统中,亥姆霍兹自由能函数包含了系统所有的热力学性质,正如下面的练习所展示的。

  • 习题:证明正则系综中系统的能量和压强的平均值可以表示为(根据量纲相等可以帮助记忆下面的公式)

\begin{equation} E = -\frac{\partial}{\partial \beta} \ln Z \end{equation} \begin{equation} p = \frac{1}{\beta} \frac{\partial}{\partial V} \ln Z \end{equation}


当然,如果系统的粒子数不是无穷大的话,计算出的热力学量应该有涨落,即标准偏差不等于零。正则系综中能量平方的平均值 \begin{equation} \langle E^2 \rangle = \sum_i w_i E_i^2 \end{equation} 可以写成 \begin{equation} \langle E^2 \rangle = \langle E \rangle^2 - \frac{\partial \langle E \rangle }{\partial\beta}. \end{equation} 其中,为了明确起见,我们把之前定义的平均能量 [math]E[/math] 写成了[math]\langle E \rangle[/math]。进而可求出能量的方差 \begin{equation} (\Delta E)^2 = \langle E^2 \rangle - \langle E \rangle^2 = - \frac{\partial \langle E \rangle }{\partial\beta} = \frac{1}{k_B \beta^2} \frac{ \partial \langle E \rangle}{\partial T} = k_B T^2 C_V. \end{equation} 其中,[math]C_V[/math] 是系统的等容热容。这就证明了等容热容一定是非负的。因为能量和等容热容都是广延量,能量的相对偏差在热力学极限(保持粒子数密度不变时让粒子数趋近于无穷大)下趋近于零: \begin{equation} \frac{\Delta E}{\langle E \rangle} = \frac{\sqrt{k_B T^2 C_V}}{\langle E \rangle} \rightarrow \frac{1}{\sqrt{N}}. \end{equation} 这个结论的数学根源就是中心极限定理的结论:即不管每一个粒子的能量如何分布,在粒子数趋近于无穷大时,系统总能量的相对偏差一定趋近于零。

巨正则系综理论

巨正则系综的概率函数

  • 巨正则系综考虑的是 [math]M[/math] 个相同的具有一定化学势、体积和温度的热力学系统。可以想象 [math]M[/math] 个相同的系统排成一排,相邻的系统之间有微弱的相互作用,使得所有的系统最终能处于同一温度。相邻的系统之间还可以交换粒子,从而每个系统的粒子数是不固定的,但根据热力学理论我们知道平衡时每个系统的化学势会相等。我们的目的是发展一套能够描述这个系综中的某一个系统的热力学性质的统计理论。大家会看到,本讲的内容与《正则系综》那一讲的内容是完全类似的。但由于巨正则系综非常重要,还是有必要完整地讨论一遍。
  • [math]M_{Ni}[/math] 表示粒子数为 [math]N[/math] 且状态为[math]i[/math] 的系统数目,而用 [math]E_{Ni}[/math] 表示对应的能量,则整个系综的总系统数、总粒子数和总能量为

\begin{equation} M = \sum_N \sum_{i(N)} M_{Ni}. \end{equation} \begin{equation} N_M = \sum_N \sum_{i(N)} M_{Ni} N. \end{equation} \begin{equation} E_M = \sum_N \sum_{i(N)} M_{Ni} E_{Ni}. \end{equation}

  • 对于一个给定的分布 [math] \{M_{Ni}\}[/math],系综的微观状态数为

\begin{equation} \Omega = \frac{M!}{\prod_N \prod_{i(N)} M_{Ni}!}. \end{equation} 从而,这个系综(作为一个很大的孤立系统)的熵为 \begin{equation} S_M = k_B \ln \Omega = k_B \ln \left(\frac{M!}{\prod_N \prod_{i(N)} M_{Ni}!}\right). \end{equation}

  • 下面我们要问:哪个分布 [math]\{M_{Ni}\}[/math] 的概率最大?那个具有最大概率的分布就对应于平衡态。根据等概率原理,肯定是微观状态数最大的分布概率最大。类比正则系综的讨论,我们知道,我们需要求 [math]\ln \Omega[/math] 的最大值。方法是引入三个拉格朗日乘子[math]\alpha[/math][math]\beta[/math][math]\gamma[/math],并将极值条件写成

\begin{equation} \frac{\partial \ln \Omega}{\partial M_{Ni}} - \alpha \frac{\partial \sum_N \sum_{j(N)} M_{Nj}}{\partial M_{Ni}} - \beta \frac{\partial\sum_N \sum_{j(N)} M_{Nj} E_{Nj}}{\partial M_{Ni}} - \gamma \frac{\partial \sum_N \sum_{j(N)} M_{Nj} N}{\partial M_{Ni}} = 0. \end{equation}

  • [math]M[/math] 趋近于无穷大时,所有的 [math]M_{Ni}[/math] 也趋近于无穷大。利用斯特林公式,可将极值条件化为

\begin{equation} \ln M_{Ni} = -\alpha - \beta E_{Ni} - \gamma N. \end{equation}

  • 对上式两边取指数,可得

\begin{equation} M_{Ni} = \exp[-\alpha - \beta E_{Ni} - \gamma N]. \end{equation} 既然 [math]M_{Ni}[/math] 是处于态 [math]Ni[/math] 的系统的个数,那么自然地,一个系统处于态 [math]Ni[/math] 的概率为 \begin{equation} w_{Ni} = \frac{M_{Ni}}{M} = \frac{\exp[- \beta E_{Ni} - \gamma N]} {\sum_N \sum_{i(N)} \exp[- \beta E_{Ni} - \gamma N]}. \end{equation}

  • 可见,常数 [math]\alpha[/math] 没有什么物理意义,在求概率的时候就被消掉了。上式中的分母叫做巨正则配分函数(grand canonical partition function),记为 [math]\Xi[/math]

\begin{equation} \boxed{\Xi = \sum_N \sum_{i(N)} \exp[- \beta E_{Ni} - \gamma N]}. \end{equation} 巨配分函数包含了体系所有的热力学性质。

巨正则系综中的热力学量

  • 根据我们对正则系综的讨论,我们可以继续认为

\begin{equation} \boxed{\beta = \frac{1}{k_B T}}. \end{equation} 这个等式的正确性可以由之后的讨论得到证实。我们不妨先相信它是正确的。下面的问题是确定系数 [math]\gamma[/math] 的物理意义。突破口还是熵和其它热力学函数。

  • 将概率函数的表达式代入熵的吉布斯公式可得

\begin{equation} S = k_B \ln \Xi + \frac{E}{T} + k_B \gamma N. \end{equation} 其中, \begin{equation} E = \sum_{N}\sum_{i(N)} w_{Ni} E_{Ni} \end{equation} 是系统能量的平均值, \begin{equation} N = \sum_{N}\sum_{i(N)} w_{Ni} N \end{equation} 是系统粒子数的平均值。 于是, \begin{equation} E - TS + k_B T \gamma N = -k_B T\ln \Xi. \end{equation} 上式右边等同于巨热力学势 \begin{equation} \boxed{\Phi = -k_B T \ln \Xi}, \end{equation} 所以,[math]\gamma[/math] 与化学式有如下联系: \begin{equation} \boxed{\gamma = -\frac{\mu}{k_B T} = -\beta \mu}. \end{equation} 于是,巨正则配分函数可以写成 \begin{equation} \boxed{\Xi = \sum_N \sum_{i(N)} \exp \left[- \frac{E_{Ni}}{k_BT} + \frac{\mu N}{k_BT} \right]}. \end{equation}


  • 容易证明,巨正则系综中系统的粒子数、能量和压强的平均值可以表示为(根据量纲相等可以帮助记忆下面的公式)

\begin{equation} \langle N \rangle = -\frac{\partial}{\partial \gamma} \ln \Xi \end{equation} \begin{equation} \langle E \rangle = -\frac{\partial}{\partial \beta} \ln \Xi \end{equation} \begin{equation} \langle p \rangle = \frac{1}{\beta} \frac{\partial}{\partial V} \ln \Xi \end{equation}

  • 当然,如果系统的粒子数不是无穷大的话,计算出的热力学量应该有涨落,即标准偏差不等于零。可以证明:巨正则系综中能量的方差为

\begin{equation} (\Delta E)^2 = \langle E^2 \rangle - \langle E \rangle^2 = - \frac{\partial \langle E \rangle }{\beta}, \end{equation} 粒子数的方差为 \begin{equation} (\Delta N)^2 = \langle N^2 \rangle - \langle N \rangle^2 = - \frac{\partial \langle N \rangle }{\gamma}. \end{equation} 于是,可以判断,在热力学极限下(即保持粒子数密度不变的条件下将粒子数增加到无穷大),能量和粒子数的相对偏差都趋近于零: \begin{equation} \frac{\Delta E}{\langle E \rangle} \rightarrow \frac{1}{\sqrt{\langle N \rangle}}. \end{equation} \begin{equation} \frac{\Delta N}{\langle N \rangle} \rightarrow \frac{1}{\sqrt{\langle N \rangle}}, \end{equation}

经典概率密度函数

  • 在得到正则配分函数之后,让我们从离散能量的情形回到连续能量的情形。为此,只要将概率函数 [math]w_i[/math] 换成概率密度函数 [math]\rho(H(q, p))[/math] 即可:

\begin{equation} \rho(H(q, p)) = \frac{e^{-\beta H(q, p)}}{\int dq dp e^{-\beta H(q, p)}}. \end{equation}

  • 这个概率密度函数显然是归一化的。这里的分母

\begin{equation} Z = \int dq dp e^{-\beta H(q, p)} \end{equation} 就是连续情形中的配分函数。这个配分函数的定义不是最终的正确形式(一个明显的问题就是 [math]Z[/math] 的量纲不为 1),但在遇到问题之前,我们不知道该对它作怎样的修正。其实,在很多的问题中,我们用不到配分函数,不妨先简单地将它看做一个归一化常数。

麦克斯韦速度与速率分布

  • 对经典粒子系统,总能量 [math]H(q, p)[/math] 是动能 [math]K(p)[/math] 和势能 [math]U(q)[/math] 的和:

\begin{equation} H(q,p) = K(p) + U(q). \end{equation} 于是,系统处于 [math]dpdq[/math] 的概率 [math]\rho(H(q, p)) dpdq[/math] 可以分解动量部分 [math]f_{p}(p)dp[/math] 和坐标部分[math]f_{q}(q)dq[/math] 的乘积: \begin{equation} \frac{e^{-\beta H(q, p)}}{Z} dpdq = \left[A e^{-\beta K(p)} dp \right] \left[B e^{-\beta U(q)} dq \right] \equiv \left[f_{p}(p)dp\right] \left[ f_{q}(q)dq\right]. \end{equation} 其中,常数 [math]A[/math][math]B[/math] 分别是动量和坐标部分的归一化因子,满足 [math]AB=1/Z[/math]。这就是说,我们可以将动量分布函数与坐标分布函数分开来研究。这其实就是关于独立随机变量的乘法原理的体现。

  • 对一般的系统,如果我们不知道体系的势能函数,是无法研究坐标分布函数的。然而,无论哪种系统(现在只讨论非相对论的经典系统;该论断在量子统计中不正确),其动量分布函数都是一样的。因为任何系统的动能都可以写成单个粒子的动能的和,所以我们可以进一步将整个系统的动量分布函数 [math]f_{p}(p)dp[/math](这里的 [math]p[/math] 指代所有[math]3N[/math] 个动量分量)分解为单个粒子的动量分布函数[math]f_{\vec{p}}(\vec{p})[/math] (这里的 [math]\vec{p}[/math] 指代某个原子的动量矢量)的乘积。其中,[math]f(\vec{p})[/math]

\begin{equation} f_{\vec{p}}(\vec{p}) dp_x dp_y dp_z = a e^{-\beta \frac{p_x^2+p_y^2+p_z^2}{2m}} dp_x dp_y dp_z = a e^{- \frac{p_x^2+p_y^2+p_z^2}{2mk_B T}} dp_x dp_y dp_z. \end{equation} 注意,此处的归一化因子 [math]a[/math] 与前面的归一化因子 [math]A[/math] 的关系为 [math]A=a^N[/math]

  • 容易计算,上述单个粒子的动量分布函数的归一化因子为

\begin{equation} a = \frac{1}{(2\pi m k_B T)^{3/2}}. \end{equation} 于是,我们可以将单粒子动量分布函数完整地写出来: \begin{equation} f_{\vec{p}}(\vec{p}) dp_x dp_y dp_z = \frac{1}{(2\pi m k_B T)^{3/2}} e^{- \frac{p_x^2+p_y^2+p_z^2}{2mk_B T}} dp_x dp_y dp_z. \end{equation}

  • 如果将随机变量从动量换成速度 [math]\vec{v}[/math],我们可以立刻写出速度分布函数:

\begin{equation} f_{\vec{v}}(\vec{v}) dv_x dv_y dv_z = \left( \frac{m}{2\pi k_B T} \right) ^{3/2} e^{- \frac{m(v_x^2+v_y^2+v_z^2)}{2k_B T}} dv_x dv_y dv_z. \end{equation} 这就是麦克斯韦于 1860 年通过分子运动论推导出的速度分布函数。

  • 我们可以继续将各个方向的速度分量的分布函数分离出来。例如,[math]x[/math] 方向的速度分布函数为

\begin{equation} f_{v_x}(v_x) dv_x = \left( \frac{m}{2\pi k_B T} \right) ^{1/2} e^{- \frac{mv_x^2}{2k_B T}} dv_x. \end{equation}

  • 如果从速度的直角坐标换到速度的球坐标([math]v[/math][math]\theta[/math][math]\phi[/math]),则有

\begin{equation} f_{\vec{v}}(\vec{v}) dv_x dv_y dv_z = \left( \frac{m}{2\pi k_B T} \right) ^{3/2} e^{- \frac{mv^2}{2k_B T}} v^2 dv \sin\theta d\theta d\phi. \end{equation} 若定义 [math]f_v(v)dv[/math] 为速率间隔 [math]dv[/math] 内的概率,则有 \begin{align} f_v(v)dv =& \left( \frac{m}{2\pi k_B T} \right) ^{3/2} e^{- \frac{mv^2}{2k_B T}} v^2 dv \int_0^{\pi} \sin\theta d\theta \int_0^{2\pi} d\phi \nonumber \\ =& 4\pi \left( \frac{m}{2\pi k_B T} \right) ^{3/2} e^{- \frac{mv^2}{2k_B T}} v^2 dv . \end{align}

  • 利用 [math]f_{v_x}(v_x)[/math] 的表达式,可以证明:

\begin{equation} \langle v_x^2 \rangle = \frac{k_B T}{m} \end{equation}

  • 利用 [math]f_{v}(v)[/math] 的表达式,可以证明:

\begin{equation} \langle v^2 \rangle = \frac{3k_B T}{m} \end{equation} \begin{equation} \langle v \rangle^2 = \frac{8k_B T}{\pi m} \end{equation}

  • 由上述结果可知,每一个平动自由度的平均能量为 [math]k_BT/2[/math]。下面我们马上会知道,这是能量均分定理的体现。

能量均分定理

  • 考虑一个经典的多粒子系统,我们用 [math]q[/math] 代表系统的 [math]s[/math] 个广义坐标,用 [math]p[/math] 代表 [math]s[/math] 个广义动量,而 [math]x_i[/math] 或者 [math]x_j[/math] 代表任意一个广义坐标或者广义动量。我们现在利用上述概率密度函数来计算经典正则系综中的一个平均值

\begin{equation} \left\langle x_i \frac{\partial H}{\partial x_j} \right\rangle = \frac{ \int dq dp e^{-\beta H(q, p)} x_i \frac{\partial H}{\partial x_j} } { \int dq dp e^{-\beta H(q, p)} }. \end{equation}

为了计算这个平均值,让我们先关注分母。注意到分母中的乘积 [math]e^{-\beta H(q, p)} \frac{\partial H}{\partial x_j}[/math] 可以写成 [math]-\frac{1}{\beta}\frac{\partial e^{-\beta H(q, p)}}{\partial x_j}[/math]。于是,分母中的被积函数可以写成 \begin{equation} -\frac{1}{\beta} x_i \frac{\partial e^{-\beta H(q, p)} }{\partial x_j} = -\frac{1}{\beta} \frac{\partial} {\partial x_j} (x_i e^{-\beta H(q, p)}) +\frac{1}{\beta} \delta_{ij} e^{-\beta H(q, p)} . \end{equation} 上式等号右边的第一项的积分等于 \begin{equation} -\frac{1}{\beta} \int \frac{dq dp}{dx_j} \int dx_j \frac{\partial} {\partial x_j} \left(x_i e^{-\beta H(q, p)}\right) = -\frac{1}{\beta} \int \frac{dq dp}{dx_j} \left( x_i e^{-\beta H(q, p)} \right)_{x_j^{min}}^{x_j^{max}} = 0 \end{equation}

  • 上式等于零的理由如下。当 [math]x_j[/math] 等于某个广义坐标时,[math]x_j^{min}[/math][math]x_j^{max}[/math] 一定对应于系统的边界处(容器壁),那里的势能无穷大,使得因子 [math]x_i e^{-\beta H(q, p)}[/math] 等于零。当 [math]x_j[/math] 等于某个广义动量时,[math]x_j^{min}[/math][math]x_j^{max}[/math] 就分别等于 [math]-\infty[/math][math]+\infty[/math],使得系统的动能为无穷大,从而依然使得因子 [math]x_i e^{-\beta H(q, p)}[/math]等于零。
  • 于是,我们要求的平均值为

\begin{equation} \left\langle x_i \frac{\partial H}{\partial x_j} \right\rangle = \frac{ \int dq dp\frac{1}{\beta} \delta_{ij} e^{-\beta H(q, p)} } { \int dq dp e^{-\beta H(q, p)} } = \frac{1}{\beta} \delta_{ij} = k_B T \delta_{ij}. \end{equation}

现在,假设系统的哈密顿量是某个广义坐标或者广义动量的二次函数,即假设由自由度 [math]x_i[/math] 贡献的哈密顿量为 \begin{equation} H_i(x_i) = a x_i^2. \end{equation} 其中,[math]a[/math] 是常数系数。对于这样的哈密顿量,显然有 \begin{equation} H_i(x_i) = \frac{1}{2} x_i \frac{\partial H_i(x_i)}{\partial x_i}. \end{equation} 于是,与自由度 [math]x_i[/math] 相关的能量平均值为 \begin{equation} \langle H_i(x_i) \rangle = \frac{1}{2} k_B T. \end{equation} 这就是能量均分定理,它说的是,哈密顿量的每个具有平方形式的自由度对总能量的平均贡献都是 [math]\frac{1}{2} k_B T[/math]

  • 将能量均分定理应用与单原子理想气体,因为每个原子只有三个平动自由度,故可得内能

\begin{equation} E = \frac{3}{2}Nk_BT. \end{equation} 于是,单原子理想气体的等容热容为 [math]C_V = \frac{3}{2}Nk_B[/math],绝热指数为 [math]\frac{3/2+1}{3/2}=5/3[/math]。该预言与实验结果是符合得很好的。

  • 再考虑双原子分子的理想气体。因为每个分子有三个平动动能项,两个转动动能项,一个振动动能项与一个振动势能项,故由能量均分定理,内能应该为

\begin{equation} E = \frac{7}{2}Nk_BT. \end{equation} 相应的等容热容为 [math]C_V = \frac{7}{2}Nk_B[/math],绝热指数为[math]\frac{7/2+1}{7/2}=9/7[/math]。然而,实验结果显示,大部分的双原子气体的绝热指数都接近于 [math]7/5[/math]。这是经典统计力学遭遇的困难之一。

经典统计力学在固体比热和黑体辐射的问题上也遭遇了巨大的困难,但这需要由量子统计解决。

位力定理

  • 很容易写出上面所定义的平均值的一个 [math]x_i=x_j=q_{\alpha}[/math] 的特例:

\begin{equation} \left\langle q_{\alpha} \frac{\partial H}{\partial q_{\alpha}} \right\rangle = k_B T. \end{equation} 如果考虑的是由 [math]N[/math] 个粒子组成的经典系统,那么,将上式对[math]\alpha[/math] 求和可得 \begin{equation} \sum_{\alpha=1}^{3N} \left\langle q_{\alpha} \frac{\partial H}{\partial q_{\alpha}} \right\rangle = 3N k_B T. \end{equation} 因为 [math]-\frac{\partial H}{\partial q_{\alpha}} = F_{\alpha}[/math] 是与坐标 [math]q_{\alpha}[/math] 对应的力,故有 \begin{equation} \sum_{\alpha=1}^{3N} \left\langle q_{\alpha} F_{\alpha} \right\rangle = -3N k_B T. \end{equation}

  • 克劳修斯于 1870 年为上式左边的量起了一个名字:Virial。常见的中文翻译有两个:一个是位力,一个是维里。我们将用第一个,因为这个量的物理意义就是位置乘以力的和的平均值,翻译成位力是非常棒的。我们将用符号 [math]W[/math] 表示这个量,即

\begin{equation} W = \sum_{\alpha=1}^{3N} \left\langle q_{\alpha} F_{\alpha} \right\rangle = -3N k_B T. \end{equation} 这个公式叫做位力定理。

现在将位力定理应用于经典理想气体,即无相互作用的经典多粒子系统。因为粒子之间没有相互作用,所有的力来自于粒子与容器的碰撞。取一个坐标系,记容器壁 [math]S[/math] 上的任意一点的坐标为 [math]\vec{x}[/math]。假设系统的压强为 [math]p[/math],则容器壁上位于 [math]\vec{x}[/math] 处的面积元 [math]d\vec{a}[/math](朝外的方向为正方向)施加给系统的力为 [math]d\vec{F} = -p d\vec{a}[/math]。于是,系统的位力可表示为 \begin{equation} W = \oint_S \vec{x} \cdot d\vec{F} = \oint_S \vec{x} \cdot (-p d\vec{a}) = -p \oint_S \vec{x} \cdot d\vec{a}. \end{equation} 记容器所在区域为 [math]\Omega[/math],包围的体积为 [math]V[/math],并应用高斯定理,可得 \begin{equation} W = -p \int_{\Omega} (\nabla \cdot \vec{x}) dv = -3p V. \end{equation} 将上式与位力定理对比可知: \begin{equation} p V = N k_B T. \end{equation} 没错,这就是理想气体的状态方程。这就证明了之前利用 [math]\beta = \frac{1}{k_B T}[/math] 所定义的 [math]T[/math] 就是绝对温度。以后我们就可以放心地将 [math]\beta = \frac{1}{k_B T}[/math] 中的[math]T[/math] 当做绝对温度了。

正则系综的应用:单原子分子理想气体

  • 虽然之前用位力定理推导了理想气体状态方程,但我们并没有系统计算单原子分子理想气体的各个热力学量。现在是时候考虑这些问题了。首先,我们明显地写出 [math]N[/math]个无相互作用的原子组成的理想气体系统的能量函数:

\begin{equation} H(q,p) = \sum_{i=1}^{N}\frac{\vec{p}_i^2}{2m} \end{equation} 对这样的系统,其配分函数是 \begin{equation} Z = \int \exp\left[-\sum_{i=1}^{N}\frac{\vec{p}_i^2}{2mk_BT} \right] \prod_{i=1}^{N} d\vec{x}_i d\vec{p}_i \end{equation}

  • 然而,正如之前就指出过的,这个配分函数的量纲都不对。要使配分函数的量纲等于 1,我们必须将上式除以一个量纲为 ([长度][math]\times[/math][动量])[math]^{3N}[/math] 的量。这样做其实就是定义一个量纲为[长度][math]\times[/math][动量]的“最小”的相空间体积 [math]\omega_0[/math],使得

\begin{equation} \int \frac{\prod_{i=1}^{N} d\vec{x}_i d\vec{p}_i}{\omega_0^{3N}} \end{equation} 等于系统中总的“相点个数”。我们期望这个最小相空间体积 [math]\omega_0[/math]应该是一个小量。

  • 你是否想到了普朗克引入的那个常数 [math]h[/math]?它的量纲与 [math]\omega_0[/math] 的量纲相同。所以,[math]\omega_0[/math] 一定正比于 [math]h[/math]。那这个比例常数是多少呢?在将由此导出的结果与其它理论结果或者实验结果对比之前,我们是无法知道答案的。暂且就假设 [math]\omega_0=h[/math] 吧。于是,配分函数变为

\begin{equation} Z = \frac{1}{h^{3N}} \int \exp\left[-\sum_{i=1}^{N}\frac{\vec{p}_i^2}{2mk_BT} \right] \prod_{i=1}^{N} d\vec{x}_i d\vec{p}_i \end{equation}

  • 容易证明,上述配分函数可以写成

\begin{equation} Z = Z_1^N. \end{equation} \begin{equation} Z_1 = \frac{V}{\lambda^3}. \end{equation} \begin{equation} \lambda =\frac{h}{\sqrt{2\pi mk_BT}}. \end{equation} 于是,体系的自由能为 \begin{equation} F = -k_BT\ln Z = -Nk_BT \ln \left( \frac{V}{\lambda^3}\right). \end{equation} 通过上述自由能,可以计算体系的压强,从而导出理想气体状态方程: \begin{equation} pV = N k_BT. \end{equation} 正如所期望的,我们推导出了理想气体应该满足的状态方程。

  • 下面再算算熵。体系的熵为:

\begin{equation} S = N k_B \ln \left( \frac{V}{\lambda^3} \right) + \frac{3}{2}Nk_B \end{equation} 虽然上述熵和自由能以及内能满足关系 [math]F=E-TS[/math],但值的注意的是熵和自由能都不是广延量。这是不可接受的。这说明我们的理论还是有不完美的地方。这个问题能由所谓的吉布斯佯谬(Gibbs paradox)更生动地展现出来。

  • 考虑由两个两个温度和数密度都相同的理想气体系统构成的孤立系统,粒子数分别为 [math]N_1[/math][math]N_2[/math]。可以证明,两个系统混合后与混合前总熵的差为

\begin{equation} \Delta S = k_B (N\ln N - N_1\ln N_1 -N_2\ln N_2 ) \approx k_B (\ln N! - \ln N_1! -\ln N_2! ). \end{equation}

  • 如果两个子系统中的气体是相同种类的气体,这个结果是很荒谬的。这就是吉布斯佯谬。同时,上述结果暗示我们,如果重新定义配分函数,使得熵的值为原来的值减去 [math]k_B \ln N![/math],也许就能消除这个佯谬。根据熵与配分函数的关系可以猜测,应该重新定义如下的配分函数

\begin{equation} Z = \frac{1}{h^{3N}N!} \int \exp\left[-\sum_{i=1}^{N}\frac{\vec{p}_i^2}{2mk_BT} \right] \prod_{i=1}^{N} d\vec{x}_i d\vec{p}_i \end{equation} 这样定义的结果就是将系统中总的相点(状态数)个数减小 [math]N![/math] 倍。重复之前的推导可得 \begin{equation} Z = \frac{Z_1^N}{N!}. \end{equation}

  • 从这个新的配分函数出发,可以证明理想气体的自由能和熵的正确表达式应该是

\begin{equation} F = -k_BT\ln Z = -Nk_BT \left[ \ln \left( \frac{v}{\lambda^3} \right) + 1 \right]. \end{equation} \begin{equation} S = N k_B \ln \left( \frac{v}{\lambda^3} \right) + \frac{5}{2}Nk_B \end{equation} 其中, \begin{equation} v = V/N. \end{equation}

  • 上述熵的公式叫做 Sackur-Tetrode 公式。正是 Sackur 在 1911 年首次假设 [math]\omega_0=h[/math],并由 Tetrode 于 1912 年最终通过将理论与实验比较确定该表达式的。由量子统计力学可以推导出该结果,但在 1911-1912 年这个时候量子力学还处于初级阶段呢。
  • 在对配分函数作出修正之后,上面计算出的自由能和熵都是广延量了。而且可以验证,吉布斯佯谬不复存在。这说明这个修正是合理的。那么,这个修正究竟代表什么意思呢?答案是:它反映了全同粒子的不可分辨性。这是量子力学中的一个结果。由于全同粒子之间不可分辨,对体系的所有粒子做一个重排并不改变系统的微观状态,故需要将体系的微观状态数在原来(没有考虑粒子的不可分辨性时)的基础上除以所有可能的重排数目,即 [math]N![/math]
  • 最后,容易计算,理想气体的化学势为

\begin{equation} \mu = -k_BT \ln \left( \frac{v}{\lambda^3}\right). \end{equation} 可见它是负的,与我们在热力学中得到结果一致。

一个简单的分子动力学模拟程序 (未修订,请不要看)

在前几章,我们回顾了经典力学、热力学及统计力学的基础知识。从本章开始,我们讨论本书的主题,即分子动力学模拟。

分子动力学模拟的定义

  • 作者曾经在一篇博文给分子动力学模拟下过一个定义:
   分子动力学模拟是一种数值计算方法,在这种方法中,我们对一个具有一定初始条件和边界条件且具有相互作用的多粒子系统的运动方程进行数值积分,得到系统在相空间中的一条离散的轨迹,并用统计力学的方法从这条相轨迹中提取出有用的物理结果。
  • 在本章余下的部分,我们会一一考察上述定义中的重要概念,如初始条件、边界条件、相互作用、运动方程、数值积分等。这里,我们首先讨论上述定义中的“统计力学的方法”。
  • 我们在第 4 章讨论过统计力学,知道在一个特定的平衡态统计系综中,一个物理量 [math]A[/math] 的统计平均值可表达为

$$ \langle A \rangle_{\rm ensemble} = \int f(p,q) A dpdq $$ 其中,[math]f(p,q)[/math] 是该系综的分布函数。这样的统计平均称为系综平均。然而,根据我们的定义,在分子动力学模拟中并没有使用“系综”(即系统的集合),而是仅对一个系统的时间演化过程进行分析。那么,上述定义中的“统计力学的方法”指的是什么呢?

  • 这里的“统计力学方法”指的是时间平均,即

$$ \langle A \rangle_{\rm time} = \lim_{t\to\infty}\frac{1}{t}\int_0^t A(t') dt' $$ 我们知道,随着时间的推移,系统将历经一条相轨迹。因为这条相轨迹永不与自身相交(第 2 章习题),如果它不代表一个周期性运动的话,那么随着时间的增加,这条相轨迹应该历经越来越多的相点。一个自然的假设是,当时间趋近于无穷大时,这条相轨迹将遍历系统的所有相点。这就是各态历经假设(ergodic hypothesis)。在此假设下,系综平均与时间平均等价。

  • 从实验的角度来说,对热力学体系的实验测量采用了时间平均而非系综平均。所以,分子动力学模拟实际上比基于系综的统计力学更加贴近实验。从计算机模拟的角度来看,除了采用时间平均的分子动力学模拟,还有直接采用系综平均的蒙特卡罗模拟,但本书不讨论蒙特卡罗模拟。

简单分子动力学模拟的基本流程

  • 根据上述定义,我们可以设想一个典型的、简单的分子动力学模拟有如下大致的计算流程:
    • 初始化。设置系统的初始条件,具体包括各个粒子的位置矢量和速度矢量。
    • 时间演化。根据系统中的粒子所满足的相互作用规律,由牛顿定律确定所有粒子的运动方程(二阶常微分方程组),并对运动方程进行数值积分,即不断地更新每个粒子的坐标和速度。最终, 我们将得到一系列离散的时刻系统在相空间中的位置,即一条离散的相轨迹。
    • 测量。用统计力学的方法分析相轨迹所蕴含的物理规律。

初始化

  • 初始化指的就是给定一个初始的相空间点,这包括各个粒子初始的坐标和速度。在分子动力学模拟中,我们需要对 [math]3N[/math][math]N[/math] 是粒子数目)个二阶常微分方程进行数值积分。因为每一个二阶常微分方程的求解都需要有两个初始条件,所以我们需要确定 [math]6N[/math] 个初始条件:[math]3N[/math] 个初始坐标分量和 [math]3N[/math] 个初始速度分量。

坐标初始化

  • 坐标的初始化指的是为系统中的每个粒子选定一个初始的位置坐标。分子动力学模拟中如何初始化位置主要取决于所要模拟的体系。例如,如果要模拟固态氩,就得让各个氩原子的位置按面心立方结构排列。如果要模拟的是液态或者气态物质,那么初始坐标的选取就可以比较随意了。重要的是,在构造的初始结构中,任何两个粒子的距离都不能太小,因为这可能导致有些粒子受到非常大的力,以至于让后面的数值积分变得非常不稳定。坐标的初始化也常被称为建模,往往需要用到一些专业的知识,例如固体物理学中的知识。


速度初始化

  • 我们知道,任何经典热力学系统在平衡时各个粒子的速度满足麦克斯韦分布。然而,作为初始条件,我们并不一定要求粒子的速度满足麦克斯韦分布。最简单的速度初始化方法是产生[math]3N[/math] 个在某个区间均匀分布的随机速度分量,再通过如下两个基本条件对速度分量进行修正。
    • 第一个条件是让系统的总动量为零。也就是说,我们不希望系统的质心在模拟的过程中跑动。分子间作用力是所谓的内力,不会改变系统的整体动量,即系统的整体动量是守恒的。只要初始的整体动量为零,在分子动力学模拟的时间演化过程中整体动量将保持为零。如果整体动量明显偏离零(相对于所用浮点数精度来说),则说明模拟出了问题。这正是判断程序是否有误的标准之一。
    • 第二个条件是系统的总动能应该与所选定的初始温度对应。我们知道,在经典统计力学中,能量均分定理成立,即粒子的哈密顿量中每一个具有平方形式的能量项的统计平均值都等于 [math]k_{\rm B} T/2[/math]。其中,[math]k_{\rm B}[/math] 是玻尔兹曼常数,[math]T[/math] 是系统的绝对温度。所以,在将质心的动量取为零之后就可以对每个粒子的速度进行一个标度变换,使得系统的初始温度与所设定的温度一致。假设我们设置的目标温度是[math]T_0[/math],那么对各个粒子的速度做如下变换即可让系统的温度从 [math]T[/math] 变成 [math]T_0[/math]

\begin{equation} \label{equation:velocity_rescale} \vec{v}_i \rightarrow \vec{v}_i'= \vec{v}_i\sqrt{\frac{T_0}{T}}. \end{equation} 容易验证,在做上式中的变换之前,如果系统的总动量已经为零,那么在做这个变换之后,系统的总动量也为零。

待写:角动量的初始化。

边界条件

  • 在我们对分子动力学模拟的定义中,除了初始条件,还提到了边界条件。边界条件对常微分方程的求解并不是必要的,但在分子动力学模拟中通常会根据所模拟的物理体系选取合适的边界条件,以期得到更合理的结果。边界条件的选取对粒子间作用力的计算也是有影响的。常用的边界条件有好几种,但我们这里只先讨论其中的一种:周期边界条件(periodic boundary conditions)。在计算机模拟中,模拟的系统尺寸一定是有限的,通常比实验中对应的体系的尺寸小很多。选取周期边界条件通常可以让模拟的体系更加接近于实际的情形,因为原本有边界的系统在应用了周期边界条件之后,“似乎”没有边界了。当然,并不能说应用了周期边界条件的系统就等价于无限大的系统,只能说周期边界条件的应用可以部分地消除边界效应,让所模拟系统的性质更加接近于无限大系统的性质。通常,在这种情况下,我们要模拟几个不同大小的系统,分析所得结果对模拟尺寸的依赖关系。
  • 在计算两个粒子,如粒子 [math]i[/math] 和粒子 [math]j[/math] 的距离时,就要考虑周期边界条件带来的影响。举个一维的例子。假设模拟在一个长度为 [math]L_x[/math] 的模拟盒子(simulation box)中进行,采用周期边界条件时,必须将该一维的盒子想象为一个圆圈。假设[math]L_x=10[/math](任意单位),第 [math]i[/math] 个粒子的坐标 [math]x_i=1[/math], 第 [math]j[/math] 个粒子的坐标 [math]x_j=8[/math],则这两个粒子的距离是多少呢?如果忽略周期边界条件,那么答案是 [math]|x_j-x_i|=7[/math],而且[math]j[/math]粒子在[math]i[/math]粒子的右边(坐标值大的一边)。但是,在采取周期边界条件时,也可认为[math]j[/math]粒子在[math]i[/math]粒子的左边,且坐标值可以平移至[math]8-10=-2[/math]。这样,[math]j[/math][math]i[/math]的距离是[math]|x_j-x_i|=3[/math],比平移[math]j[/math]粒子之前两个粒子之间的距离要小。在我们的模拟中,总是采用最小镜像约定(minimum image convention):在计算两个粒子的距离时,总是取最小的可能值。定义

\begin{equation} x_j-x_i \equiv x_{ij} \end{equation} 则这个约定等价于如下规则:

    • 如果 [math]x_{ij}\lt -L_x/2[/math],则将[math]x_{ij}[/math]换为[math]x_{ij}+L_x[/math]
    • 如果 [math]x_{ij}\gt +L_x/2[/math],则将[math]x_{ij}[/math]换为[math]x_{ij}-L_x[/math]

最终效果就是让变换后的[math]x_{ij}[/math]的绝对值小于[math]L_x/2[/math]

  • 很容易将上述讨论推广到二维和三维的情形。例如,在二维的情形中,就要将一个周期的模拟盒子想象为一个环面(torus),就像一个甜甜圈或一个充了气的轮胎。在三维的情形中,就要将一个周期的模拟盒子想象为一个三维环面,而最小镜像约定可以表达为:
    • 如果 [math]x_{ij}\lt -L_x/2[/math],则将[math]x_{ij}[/math]换为[math]x_{ij}+L_x[/math]
    • 如果 [math]x_{ij}\gt +L_x/2[/math],则将[math]x_{ij}[/math]换为[math]x_{ij}-L_x[/math]
    • 如果 [math]y_{ij}\lt -L_y/2[/math],则将[math]y_{ij}[/math]换为[math]x_{ij}+L_x[/math]
    • 如果 [math]y_{ij}\gt +L_y/2[/math],则将[math]y_{ij}[/math]换为[math]y_{ij}-L_y[/math]
    • 如果 [math]z_{ij}\lt -L_z/2[/math],则将[math]z_{ij}[/math]换为[math]z_{ij}+L_z[/math]
    • 如果 [math]z_{ij}\gt +L_z/2[/math],则将[math]z_{ij}[/math]换为[math]z_{ij}-L_z[/math]

这里,我们假设了三维模拟盒子中3个共点的边的长度分别为[math]L_x[/math][math]L_y[/math][math]L_z[/math],且两两相互垂直(所谓的正交模拟盒子)。如果有任意两个共点的边不是相互垂直的,情况就要复杂一些。本章仅讨论正交盒子的情形,以后再讨论非正交盒子的情形。


相互作用

两体势

宏观物质的性质在很大程度上是由微观粒子之间的相互作用力决定的。所以,对粒子间相互作用力的计算在分子动力学模拟中是至关重要的。粒子间有何种相互作用不是分子动力学模拟本身所能回答的;它本质上是一个量子力学的问题。在经典分子动力学模拟中,粒子间的相互作用力常常由一个或多个经验势函数描述。经验势函数能够在某种程度上反映出某些物质的某些性质。现在已发展出很多这样的势函数。这里,我们只介绍本章将要用到的一种势函数——Lennard-Jones势。考虑系统中的任意粒子对[math]i[/math][math]j[/math],它们之间的相互作用势能可以写为 \begin{equation} \label{equation:Uij} U_{ij}(r_{ij})=4\epsilon \left( \frac{\sigma^{12}}{r_{ij}^{12}}-\frac{\sigma^{6}}{r_{ij}^{6}} \right). \end{equation} 其中,[math]\epsilon[/math][math]\sigma[/math]是势函数中的参量,分别具有能量和长度的量纲;[math]r_{ij} = |\vec{r}_j - \vec{r}_i|[/math]是两个粒子间的距离。

Lennard-Jones势比较适合描述惰性元素组成的物质。它是最早提出的两体势函数之一。所谓两体势,指的是两个粒子[math]i[/math][math]j[/math]之间的相互作用势仅依赖于它们之间的距离[math]r_{ij}[/math],不依赖于系统中其他粒子的存在与否及具体位置。对于这样的势函数,我们可以将整个系统的势能[math]U[/math]写为 \begin{equation} \label{equation:U} U=\sum_{i=1}^N U_i; \end{equation} \begin{equation} \label{equation:U_i} U_i= \frac{1}{2} \sum_{j \neq i} U_{ij}(r_{ij}). \end{equation} 将以上两式合起来,可以写成 \begin{equation} U=\frac{1}{2}\sum_{i=1}^N \sum_{j \neq i} U_{ij}(r_{ij}) \end{equation} 上面的[math]U_i[/math]可以称为粒子[math]i[/math]的势能。也可以将总势能写为如下形式: \begin{equation} \label{equation:U(i<j)} U=\sum_{i=1}^N \sum_{j > i} U_{ij}(r_{ij}) \end{equation}

练习:两体势中粒子受力的推导

(1) 对于两体势,证明粒子 [math]i[/math]所受总的力为 \begin{equation} \vec{F}_{i} = \sum_{j \neq i} \vec{F}_{ij} \end{equation} \begin{equation} \vec{F}_{ij} = \frac{\partial U_{ij}(r_{ij})}{\partial r_{ij}} \frac{\vec{r}_{ij} }{r_{ij}} \end{equation} 其中,我们定义了一个表示粒子间相对位置的符号 \begin{equation} \vec{r}_{ij} \equiv \vec{r}_j - \vec{r}_i \end{equation}

(2)证明牛顿第三定律成立: \begin{equation} \vec{F}_{ij} = - \vec{F}_{ji}. \end{equation}

(3) 证明对于Lennard-Jones势,有 \begin{equation} \label{equation:Fij} \vec{F}_{ij} = \left( \frac{24 \epsilon \sigma^6}{r_{ij}^8} - \frac{48 \epsilon \sigma^{12}}{r_{ij}^{14}} \right) \vec{r}_{ij}. \end{equation}

近邻列表

Verlet列表

Cell list算法

运动方程的数值积分

The aim of time evolution is to find the phase trajectory \begin{equation} \{ \vec{r}_i(t_1), ~\vec{v}_{i}(t_1)\}_{i=1}^N,~ \{ \vec{r}_i(t_2), ~\vec{v}_{i}(t_2)\}_{i=1}^N,~ \cdots \end{equation} starting from the initial phase point \begin{equation} \{ \vec{r}_i(t_0), ~\vec{v}_{i}(t_0)\}_{i=1}^N. \end{equation} The time interval between two time points $\Delta t=t_1-t_0=t_2-t_1=\cdots$ is called the time step.


In the $NVE$ ensemble, the dynamics of the system is Hamiltonian and the equations of motion can be derived from Hamilton's equations. Because these equations of motion have the time-reversal symmetry, a good numerical integrating method (an integrator) should preserve this symmetry.

One of the most widely used integrators which has the property of time-reversibility is the so-called velocity-Verlet method . This integrator is also symplectic. These two properties make the velocity-Verlet integrator very stable for long-time simulations. Here are the velocity and position updating equations in the velocity-Verlet method: \begin{equation} \label{equation:velocity-Verlet-velocity} \vec{v}_i(t_{m+1}) \approx \vec{v}_i(t_{m}) + \frac{\vec{F}_i(t_m)+\vec{F}_i(t_{m+1})}{2m_i}\Delta t; \end{equation} \begin{equation} \label{equation:velocity-Verlet-position} \vec{r}_i(t_{m+1}) \approx \vec{r}_i(t_{m}) + \vec{v}_i(t_m) \Delta t + \frac{1}{2} \frac{\vec{F}_i(t_m)}{m_i} (\Delta t)^2, \end{equation} where $m_i$ is the mass of particle $i$.

The above velocity-Verlet integrator can be derived by finite-difference method (Taylor series expansion), but a more general method, which can be generalized to more sophisticated situations, is the classical time-evolution operator approach, or the Liouville operator approach. In this approach, the time-evolution of a classical system by one step can be formally expressed as \begin{equation} \left( \begin{array}{c} \vec{r}_i(t+\Delta t) \\ \vec{p}_i(t+\Delta t) \end{array} \right) = e^{iL\Delta t} \left( \begin{array}{c} \vec{r}_i(t) \\ \vec{p}_i(t) \end{array} \right), \end{equation} where $\vec{p}_i$ is the momentum of particle $i$ and $e^{iL\Delta t}$ is called the classical evolution operator, which is the classical counterpart of the quantum evolution operator. The operator $iL$ in the exponent of the evolution operator is called the Liouville operator and is defined by \begin{equation} iL (\text{anything}) = \{\text{anything}, H\} \equiv \sum_{i=1}^N \left( \frac{\partial H}{\partial \vec{p}_i} \cdot \frac{\partial }{\partial \vec{r}_i} - \frac{\partial H}{\partial \vec{r}_i} \cdot \frac{\partial }{\partial \vec{p}_i} \right) (\text{anything}). \end{equation} Here, $H$ is the Hamiltonian of the system. Because \begin{equation} \frac{\partial H}{\partial \vec{p}_i} = \frac{\vec{p}_i}{m_i} ~ \text{and} ~ -\frac{\partial H}{\partial \vec{r}_i} = \vec{F}_i, \end{equation} we have \begin{equation} iL = iL_1 + iL_2, \end{equation} \begin{equation} iL_1 = \sum_{i=1}^N \frac{\vec{p}_i}{m_i} \cdot \frac{\partial }{\partial \vec{r}_i}, \end{equation} \begin{equation} iL_2 = \sum_{i=1}^N \vec{F}_i \cdot \frac{\partial }{\partial \vec{p}_i}. \end{equation} Here, we have divided the Liouville operator into two parts. In general, $iL_1$ and $iL_2$ do not commute, and therefore $e^{iL \Delta t} \neq e^{iL_1 \Delta t}e^{iL_2 \Delta t}$. However, there is an important theorem called the Trotter theorem, which can be used to derive the following approximation: \begin{equation} e^{iL \Delta t} \approx e^{iL_2 \Delta t/2} e^{iL_1 \Delta t}e^{iL_2 \Delta t/2}. \end{equation} Now, we can express the one-step integration as \begin{equation} \left( \begin{array}{c} \vec{r}_i(t+\Delta t) \\ \vec{p}_i(t+\Delta t) \end{array} \right) \approx e^{iL_2 \Delta t/2} e^{iL_1 \Delta t}e^{iL_2 \Delta t/2} \left( \begin{array}{c} \vec{r}_i(t) \\ \vec{p}_i(t) \end{array} \right). \end{equation} To make further derivations, we note that for an arbitrary constant $c$, we have \begin{equation} e^{c \frac{\partial}{\partial x}} x = x+c. \end{equation} Applying this identity to the right most operator in the above equation, we have \begin{equation} \left( \begin{array}{c} \vec{r}_i(t+\Delta t) \\ \vec{p}_i(t+\Delta t) \end{array} \right) \approx e^{iL_2 \Delta t/2} e^{iL_1 \Delta t} \left( \begin{array}{c} \vec{r}_i(t) \\ \vec{p}_i(t) + \frac{\Delta t}{2} \vec{F}_i(t) \end{array} \right). \end{equation} Then, applying the operator $e^{iL_1 \Delta t}$, we have \begin{equation} \left( \begin{array}{c} \vec{r}_i(t+\Delta t) \\ \vec{p}_i(t+\Delta t) \end{array} \right) \approx e^{iL_2 \Delta t/2} \left( \begin{array}{c} \vec{r}_i(t) + \Delta t \frac{\vec{p}_i(t) + \frac{\Delta t}{2} \vec{F}_i(t)}{m_i} \\ \vec{p}_i(t) + \frac{\Delta t}{2} \vec{F}_i(t) \end{array} \right). \end{equation} Last, applying the remaining operator $e^{iL_2 \Delta t/2}$, we have \begin{equation} \left( \begin{array}{c} \vec{r}_i(t+\Delta t) \\ \vec{p}_i(t+\Delta t) \end{array} \right) \approx \left( \begin{array}{c} \vec{r}_i(t) + \Delta t \frac{\vec{p}_i(t) + \frac{\Delta t}{2} \vec{F}_i(t)}{m_i} \\ \vec{p}_i(t) + \frac{\Delta t}{2} \vec{F}_i(t) + \frac{\Delta t}{2} \vec{F}_i(t+\Delta t) \end{array} \right). \end{equation} It is clear that this equation is equivalent to Eqs. (\ref{equation:velocity-Verlet-velocity}) and (\ref{equation:velocity-Verlet-position}).


We can see that in the velocity-Verlet integrator, the position updating can be done in one step, but the velocity updating can only be done by two steps, one before force updating and the other after it.


我们知道在经典力学中,粒子的运动方程可以用牛顿第二定律表达。例如,对于第[math]i[/math]个粒子,其运动方程为 \begin{equation} m_i \frac{d\vec{r}_i^2}{dt^2} = \vec{F}_i \end{equation} 这是一个二阶常微分方程,我们可以把它改写为两个一阶常微分方程: \begin{equation} \label{equation:velocity_definition} \frac{d\vec{r}_i}{dt} = \vec{v}_i \end{equation} \begin{equation} \frac{d\vec{v}_i}{dt} = \frac{\vec{F}_i}{m_i} \end{equation}

对运动方程进行数值积分的目的就是在给定的初始条件下找到各个粒子在一系列离散的时间点的坐标和速度值。我们假设每两个离散的时间点之间的间隔是固定的,记为[math]\Delta t[/math],称为时间步长(time step)。在分子动力学模拟中使用的数值积分方法有很多种,本章只介绍所谓的“速度-Verlet”积分方法(Verlet 是发展分子动力学模拟方法的先驱之一),而且不去追究其推导过程。在该方法中,第[math]i[/math]个粒子在时刻[math]t+\Delta t[/math]的速度[math]\vec{v}_i(t+\Delta t)[/math]和位置[math]\vec{r}_i(t+\Delta t)[/math]分别由以下两式给出: \begin{equation} \vec{v}_i(t+\Delta t)=\vec{v}_i(t)+\frac{1}{2}\frac{\vec{F}_i(t)+\vec{F}_i(t+\Delta t)}{m_i}\Delta t \end{equation} \begin{equation} \vec{r}_i(t+\Delta t) =\vec{r}_i(t) +\vec{v}_i(t)\Delta t +\frac{1}{2}\frac{\vec{F}_i(t)}{m_i}(\Delta t)^2 \end{equation} 由以上两式可以看出,[math]t+\Delta t[/math]时刻的坐标仅依赖于[math]t[/math]时刻的坐标、速度和力,但[math]t+\Delta t[/math]时刻的速度依赖于[math]t[/math]时刻的速度、力及[math]t+\Delta t[/math]时刻的力。所以,从算法的角度来说,以上两式应该对应如下的计算流程: \begin{equation} \label{equation:verlet1v} \vec{v}_i(t) \rightarrow \vec{v}_i(t+\Delta t/2)=\vec{v}_i(t)+\frac{1}{2}\frac{\vec{F}_i(t)}{m_i}\Delta t \end{equation} \begin{equation} \label{equation:verlet1r} \vec{r}_i(t)\rightarrow \vec{r}_i(t+\Delta t) =\vec{r}_i(t) +\vec{v}_i(t+\Delta t/2)\Delta t \end{equation} \begin{equation} \label{equation:force_update} \vec{F}_i(t)\rightarrow \vec{F}_i(t+\Delta t) \end{equation} \begin{equation} \label{equation:verlet2v} \vec{v}_i(t+\Delta t/2) \rightarrow \vec{v}_i(t+\Delta t)=\vec{v}_i(t+\Delta t/2)+\frac{1}{2}\frac{\vec{F}_i(t+\Delta t)}{m_i}\Delta t \end{equation} 注意,我们引入了一个中间的速度变量[math]\vec{v}_i(t+\Delta t/2)[/math]。完成上述计算之后,粒子的坐标、速度、和力都从[math]t[/math]时刻的更新为[math]t+\Delta t[/math]时刻的。这就是一个时间步的计算。反复执行这样的计算流程,系统的微观状态就会不断地随时间变化,从而得到一条相空间的轨迹。系统所有的宏观性质都包含在相轨迹中。

C++ 开发一个简单的分子动力学模拟程序

  • 本节给出一个简单的用 C++ 开发的分子动力学模拟程序。我们选择用 C++ 语言开发程序,因为这是作者最熟悉的编程语言。因为分子动力学模拟一般较为耗时,所以用高效的编译型语言开发较为合适。除了 C++ 语言,C 语言和 Fortran 语言也很高效。如果需要的话,读者可使用其中的任何一门语言编写分子动力学模拟程序。另外,在完成分子动力学模拟的计算后,我们往往需要对得到的数据进行可视化,此时使用一门解释性语言更为方便。我们选择使用流行的 Python 语言进行数据的后处理和可视化。如果读者有 MATLAB 软件的话,使用它进行数据的后处理和可视化也不错。

程序的结构

程序的编译与运行

我们用 make 程序从源文件构建出可执行文件。为此,我们准备了需要的Makefile文件。如果要在Windows系统使用的话,可以从 该网站 下载一个GNU Make For Windows。可用以下命令编译:

make

程序通过编译后,会产生一个名为ljmd的二进制文件。

程序中使用的单位制

  • 我们的分子动力学模拟程序只涉及经典力学和热力学,故只需要用到4个基本物理量的单位。我们选择如下4个基本单位来确定各个物理量的数值:
  • 能量:电子伏特(记号为 eV),约为[math]1.6\times 10^{-19}[/math] J。
  • 长度:埃(angstrom,记号为A),即[math]10^{-10}[/math] m。
  • 质量:原子质量单位(atomic mass unit,记号为amu),约为[math]1.66 \times 10^{-27}[/math] kg。
  • 温度:开尔文(记号为 K)。

用这样的基本单位,可使程序中大部分物理量的数值都接近[math]1[/math]。我们称这样的单位为该程序的“自然单位”。

  • 从以上基本单位可以推导出程序中其他相关物理量的单位:
  • 力。因为力乘以距离等于功(能量),故力的单位是能量单位除以长度单位,即 eV A[math]^{-1}[/math]
  • 速度。因为动能正比于质量乘以速度的平方,故速度的单位是能量单位除以质量单位再开根号,即eV[math]^{1/2}[/math] amu[math]^{-1/2}[/math]
  • 时间。因为长度等于速度乘以时间,故时间的单位是长度单位除以速度单位,即A amu[math]^{1/2}[/math] eV[math]^{-1/2}[/math],约为[math]1.018051 \times 10^{1}[/math] fs(飞秒,即[math]10^{-15}[/math] s)。在程序的初始化阶段,我们需要设置一个积分的时间步长。一般来说,我们习惯用飞秒为单位设置积分步长。所以,我们需要随后将积分步长的单位从飞秒转换为该程序的自然单位,这需要除以如下常数:
 const double TIME_UNIT_CONVERSION = 1.018051e+1;

该常数在头文件common.h中定义。

  • 玻尔兹曼常数[math]k_{\rm B}[/math]。这是一个很重要的常数,它在国际单位制中约为[math]1.38\times 10^{-23}[/math] J K[math]^{-1}[/math],对应于程序自然单位制的 [math]8.617343 \times 10^{-5}[/math] eV K[math]^{-1}[/math]。所以,我们在头文件common.h 还定义了一个常数:
 const double K_B = 8.617343e-5;


程序正确性的测试

能量守恒的测试

本章程序不是一次性开发出来的,而是一点一点地写出来的,而且作者在编写的过程中做过很多测试,也犯过很多错误。这个开发过程是很难在书中体现出来的。我们首先测试该程序是否能通过能量守恒的测试。

程序在产出阶段每隔100步输出系统的总动能[math]K(t)[/math]和总势能[math]U(t)[/math],它们都是时间[math]t[/math](从产出阶段开始计时)的函数。对于大小有限的体系,它们都是随时间[math]t[/math]涨落的。然而,根据能量守恒定律,系统动能和势能的和,即总能量[math]E(t)=K(t)+V(t)[/math],应该是不随时间变化的。当然,我们的模拟中使用了具有一定误差的数值积分方法,故总能量也会有一定大小的涨落。这个涨落主要与积分的时间步长有关系。一般来说,积分的时间步长越大,总能量的涨落越大。

图xxx(a)给出了系统的总动能、总势能和总能量在产出过程中随时间变化的情况。可以看出动能是正的,势能是负的,涨落相对较大;总能是负的,但在该图中看不出有涨落。图xxx(b)给出了[math](E(t)-\langle E\rangle)/|\langle E\rangle|[/math],即总能的相对涨落值。可见总能量确实也有涨落,相对涨落值在[math]10^{-5}[/math]量级。对于很小的体系来说,这是一个合理的值。

练习:动量守恒和速度分布的测试

对能量守恒的检验只是判断一个分子动力学模拟程序是否正确的方法之一。一个分子动力学模拟程序通过了能量守恒的检验,不代表就没有错误了。请读者修改本章的程序,以某一个时间间隔输出各个粒子的速度,然后检验如下两点: (1)体系的动量受否守恒? (2)粒子的速度是否满足麦克斯韦速度分布规律?


统计系综与积分算法(未修订,请不要看)

NVT 系综与相关积分算法

Berendsen

Using the Berendsen thermostat, the integration algorithm in the $NVT$ ensemble only requires an extra scaling of all the velocity components. For the $NPT$ ensemble, the Berendsen barostat requires an extra scaling of positions and box lengths. The Berendsen thermostat and barosat are very suitable for equilibrating the system to a target temperature and pressure.


The velocities are scaled in the Berendsen thermostat in the following way: \begin{equation} \vec{v}_i^{\text{scaled}} = \vec{v}_i \sqrt{1 + \alpha_T \left(\frac{T_0}{T} - 1\right)}. \end{equation} Here, $\alpha_T$ is a dimensionless parameter, $T_0$ is the target temperature, and $T$ is the instant temperature calculated from the current velocities $\{ \vec{v}_i \}$. The parameter $\alpha_T$ should be positive and not larger than 1. When $\alpha_T=1$, the above formula reduces to the simple velocity-scaling formula: \begin{equation} \vec{v}_i^{\text{scaled}} = \vec{v}_i \sqrt{\frac{T_0}{T}}. \end{equation} A smaller $\alpha_T$ represents a weaker coupling between the system and the thermostat. Practically, any value of $\alpha_T$ in the range of $0.001 \sim 1$ can be used.

Nosé-Hoover

Nosé-Hoover chain

The Nos\'e-Hoover chain method is more suitable for calculating equilibrium properties in a specific ensemble. In the current version of GPUMD, only the Nos\'e-Hoover chain thermostat is implemented. We hope to implement the Nos\'e-Hoover chain barostat in a future version.

The equations of motion in the Nos\'{e}-Hoover chain method are \begin{equation} \frac{d}{dt} \vec{r}_i = \frac{\vec{p}_i}{m_i}, \end{equation} \begin{equation} \frac{d}{dt} \vec{p}_i = \vec{F}_i - \frac{\pi_0}{Q_0} \vec{p}_i, \end{equation} \begin{equation} \frac{d}{dt} \eta_k = \frac{\pi_k}{Q_k} ~(k = 0, 1, \cdots, M-1), \end{equation} \begin{equation} \frac{d}{dt} \pi_0 = 2\left( \sum_i \frac{\vec{p}_i^2}{2m_i} - dN\frac{k_BT}{2} \right) - \frac{\pi_1}{Q_1} \pi_0, \end{equation} \begin{equation} \frac{d}{dt} \pi_k = 2\left( \frac{\pi_{k-1}^2}{2Q_{k-1}} - \frac{k_BT}{2} \right) - \frac{\pi_{k+1}}{Q_{k+1}} \pi_{k} ~(k = 1, 2, \cdots, M-2), \end{equation} \begin{equation} \frac{d}{dt} \pi_{M-1} = 2\left( \frac{\pi_{M-2}^2}{2Q_{M-2}} - \frac{k_BT}{2} \right). \end{equation} The optimal choice for the thermostat masses is \begin{equation} Q_0 = dNk_BT\tau^2, \end{equation} \begin{equation} Q_k = k_BT\tau^2 ~(k = 1, 2, \cdots, M-1), \end{equation} where $\tau$ is a time parameter, whose value is usually chosen by try and error in practice. A good choice is $\tau = 100 \Delta t$, where $\Delta t$ is the time step for integration.

An integration scheme for the $NVT$ ensemble using the Nos\'{e}-Hoover chain can also be formulated using the approach of the time-evolution operator. The total Liouville operator for the equations of motion in the Nos\'{e}-Hoover chain method is \begin{equation} iL = iL_1 + iL_2 + iL_T, \end{equation} \begin{equation} iL_1 = \sum_{i=1}^N \frac{\vec{p}_i}{m_i} \cdot \frac{\partial }{\partial \vec{r}_i}, \end{equation} \begin{equation} iL_2 = \sum_{i=1}^N \vec{F}_i \cdot \frac{\partial }{\partial \vec{p}_i}. \end{equation} \begin{equation} iL_T = \sum_{k=0}^{M-1} \frac{\pi_k}{Q_k}\frac{\partial}{\partial \eta_k} + \sum_{k=0}^{M-2} \left(G_k - \frac{\pi_{k+1}}{Q_{k+1}}\pi_k\right)

    \frac{\partial}{\partial \pi_k}                               +

G_{M-1} \frac{\partial}{\partial \pi_{M-1}} - \sum_{i=0}^{N-1} \frac{\pi_0}{Q_0} \vec{p}_i \cdot \frac{\partial}{\partial \vec{p}_i}. \end{equation} That is, the Liouville operator for the $NVT$ ensemble contains an extra term $iL_T$ related to the thermostat variables, which is absent from that for the $NVE$ ensemble.

The total time-evolution operator $e^{iL\Delta t} $ for one step can be factorized using the Trotter theorem as in the case of the $NVE$ ensemble: \begin{equation} e^{iL} \approx e^{iL_T\Delta t/2} e^{iL_2\Delta t/2} e^{iL_1\Delta t} e^{iL_2\Delta t/2} e^{iL_T\Delta t/2}. \end{equation} Comparing this with the factorization in the $NVE$ ensemble, we see that we only need to apply the operator $e^{iL_T\Delta t/2}$ before and after applying the usual velocity-Verlet integrator in the $NVE$ ensemble.

The operator $e^{iL_T\Delta t/2}$ can be further factorized into some elementary factors using the Trotter theorem. First, we define the following decomposition of the operator $iL_T$: \begin{equation} iL_T = iL_{T1} + iL_{T2} + iL_{T3}, \end{equation} \begin{equation} iL_{T1} = \sum_{k=0}^{M-1} \frac{\pi_k}{Q_k}\frac{\partial}{\partial \eta_k}, \end{equation} \begin{equation} iL_{T2} = \sum_{k=0}^{M-2} \left(G_k - \frac{\pi_{k+1}}{Q_{k+1}}\pi_k\right)

    \frac{\partial}{\partial \pi_k}                               +

G_{M-1} \frac{\partial}{\partial \pi_{M-1}}, \end{equation} \begin{equation} iL_{T3} = -\sum_{i=0}^{N-1} \frac{\pi_0}{Q_0} \vec{p}_i \cdot \frac{\partial}{\partial \vec{p}_i}. \end{equation} We can then make the following factorization: \begin{equation} e^{iL_T\Delta t/2} \approx e^{iL_{T2}\Delta t/4} e^{iL_{T3}\Delta t/2} e^{iL_{T1}\Delta t/2} e^{iL_{T2}\Delta t/4}. \end{equation} There are still a few terms in $iL_{T2}$ and we need to factorize $e^{iL_{T2}\Delta t/4}$ further. We can factorize the $e^{iL_{T2}\Delta t/4}$ term on the right of the above equation as \begin{equation} e^{iL_{T2}\Delta t/4} \approx \prod_{k=0}^{M-2} \left( e^{-\frac{\Delta t}{8} \frac{\pi_{k+1}}{Q_{k+1}}\pi_k \frac{\partial}{\partial \pi_k} } e^{\frac{\Delta t}{4} G_k \frac{\partial}{\partial \pi_k} } e^{-\frac{\Delta t}{8} \frac{\pi_{k+1}}{Q_{k+1}}\pi_k \frac{\partial}{\partial \pi_k} } \right) e^{\frac{\Delta t}{4} G_{M-1} \frac{\partial}{\partial \pi_{M-1}} } \end{equation} and correspondingly factorize that on the left as \begin{equation} e^{iL_{T2}\Delta t/4} \approx e^{\frac{\Delta t}{4} G_{M-1} \frac{\partial}{\partial \pi_{M-1}} } \prod_{k=M-2}^{0} \left( e^{-\frac{\Delta t}{8} \frac{\pi_{k+1}}{Q_{k+1}}\pi_k \frac{\partial}{\partial \pi_k} } e^{\frac{\Delta t}{4} G_k \frac{\partial}{\partial \pi_k} } e^{-\frac{\Delta t}{8} \frac{\pi_{k+1}}{Q_{k+1}}\pi_k \frac{\partial}{\partial \pi_k} } \right). \end{equation}

It can be shown that the effect of the operator $e^{cx\frac{\partial}{\partial x}}$ on $x$ is to scale it by a factor of $e^c$: \begin{equation} e^{cx\frac{\partial}{\partial x}} x = e^c x. \end{equation} Therefore, the effect of the operator $e^{iL_{T3}\Delta t/2}$ is to scale the momenta of all the particles in the system by a uniform factor $e^{-(\pi_0/Q_0)\Delta t/2}$. Although this operator appears in the factorization of $e^{iL_{T}\Delta t/2}$, its does not affect the thermostat variables. Therefore, when applying the operator $e^{iL_{T}\Delta t/2}$, we only need to update the variables related to the thermostats and save this factor for later use when we update the variables for the particles. In this way, the update for the thermostat variables and that for the particle variables are separated.

BDP

Langevin

NPT系综与相关积分算法

位力应力的详细推导

Berendsen 控压算法

In the Berendsen barostat algorithm, the particle positions and box length in a given direction are scaled if periodic boundary conditions are applied to that direction. The scaling of the positions reads \begin{equation} \vec{r}_i^{\text{scaled}} = \vec{r}_i \left[ 1-\alpha_p (\vec{p}_0 - \vec{p}) \right]. \end{equation} Here, $\alpha_p$ is a parameter and $\vec{p}_0$ ($\vec{p}$) is the target (instant) pressure in the three directions. The parameter $\alpha_p$ is not dimensionless, and it requires some try-and-error to find a good value of it for a given system. A harder/softer system requires a smaller/larger value of $\alpha_p$. In the unit system adopted by GPUMD, it is recommended that $\alpha_p = 10^{-4} \sim 10^{-2}$. Only directions with periodic boundary conditions will be affected by the barostat.

MTTK 控压算法

In this ensemble, we introduce a box matrix $\mathbf{h}$: \begin{equation} \mathbf{h} = \left( \begin{array}{ccc} a_{x} & b_{x} & c_{x} \\ a_{y} & b_{y} & c_{y} \\ a_{z} & b_{z} & c_{z} \\ \end{array} \right) \end{equation} Generally, there are 6 nonzero matrix elements in the box matrix, but in the current version of GPUMD, we only consider the special case where the box matrix is diagonal. This corresponds to orthogonal box. We may consider the general triclinic box in the future. However, the formulations below are presented in a general form.

\textbf{Important}: we use bold face to represent (second rank) tensors and italic bold face to represent vectors. The dot product of a tensor and a vector is a vector and the dot product of two tensors is a tensor. The unit tensor is denoted as $\mathbf{I}$. The tensors are also regarded as matrices and various matrix operators (such as determinant, trace, transpose, and inverse) are used.

List of dynamic variables involved: \begin{itemize} \item particle positions: $\vec{r}_i$ \item particle momenta: $\vec{p}_i$ \item box matrix: $\mathbf{h}$ \item box momenta: $\mathbf{p}_g$ \item particle-thermostat positions: $\eta_k$ \item particle-thermostat momenta: $\pi_k$ \item box-thermostat positions: $\xi_k$ \item box-thermostat momenta: $\beta_k$ \end{itemize}

The equations of motion for the particles are:

\begin{equation} \frac{d}{dt} \vec{r}_i = \frac{\vec{p}_i}{m_i} + \frac{\mathbf{p}_g}{W_g} \cdot \vec{r}_i; \end{equation}

\begin{equation} \frac{d}{dt} \vec{p}_i = \vec{F}_i - \frac{\pi_0}{Q_0} \vec{p}_i - \frac{\mathbf{p}_g}{W_g} \cdot \vec{p}_i - \frac{1}{N_f}\frac{\textmd{Tr}[\mathbf{p}_g]}{W_g} \vec{p}_i; \end{equation}

The equations of motion for the box are: \begin{equation} \frac{d\mathbf{h}}{dt} = \frac{\mathbf{p}_g}{W_g} \cdot \mathbf{h}; \end{equation} \begin{equation} \frac{d\mathbf{p}_g}{dt} = \det [\mathbf{h}] (\boldsymbol{\sigma}-\boldsymbol{\sigma}^{\rm ext}) + \frac{1}{N_f} \sum_i \frac{\vec{p}_i^2}{2m_i}\mathbf{I} -\frac{\pi'_0}{Q'_0} \mathbf{p}_g; \end{equation}

The equations of motion for the thermostat chain coupling to the particles are:

\begin{equation} \frac{d}{dt} \eta_k = \frac{\pi_k}{Q_k} ~(k = 0, 1, \cdots, M-1), \end{equation}

\begin{equation} \frac{d}{dt} \pi_0 = 2\left( \sum_i \frac{\vec{p}_i^2}{2m_i} - dN\frac{k_BT}{2} \right) - \frac{\pi_1}{Q_1} \pi_0, \end{equation}

\begin{equation} \frac{d}{dt} \pi_k = 2\left( \frac{\pi_{k-1}^2}{2Q_{k-1}} - \frac{k_BT}{2} \right) - \frac{\pi_{k+1}}{Q_{k+1}} \pi_{k} ~(k = 1, 2, \cdots, M-2), \end{equation}

\begin{equation} \frac{d}{dt} \pi_{M-1} = 2\left( \frac{\pi_{M-2}^2}{2Q_{M-2}} - \frac{k_BT}{2} \right). \end{equation}

The equations of motion for the thermostat chain coupling to the box are:

\begin{equation} \frac{d}{dt} \xi_k = \frac{\beta_k}{Q'_k} ~(k = 0, 1, \cdots, M-1), \end{equation}

\begin{equation} \frac{d}{dt} \beta_0 = 2\left(

\frac{\textmd{Tr} [\mathbf{p}_g^{T} \mathbf{p}_g] }{2W_g} - (d^2-d)\frac{k_BT}{2}

\right) - \frac{\beta_1}{Q'_1} \beta_0, \end{equation}

\begin{equation} \frac{d}{dt} \beta_k = 2\left( \frac{\beta_{k-1}^2}{2Q'_{k-1}} - \frac{k_BT}{2} \right) - \frac{\beta_{k+1}}{Q'_{k+1}} \beta_{k} ~(k = 1, 2, \cdots, M-2), \end{equation}

\begin{equation} \frac{d}{dt} \beta_{M-1} = 2\left( \frac{\beta_{M-2}^2}{2Q'_{M-2}} - \frac{k_BT}{2} \right). \end{equation}


The optimal choice for the box mass is \begin{equation} W_g = (dN+1)k_BT\tau_b^2, \end{equation} where $\tau_b$ is a time parameter, whose value is usually chosen by try and error in practice. A good choice is $\tau_b = 1000 \Delta t$, where $\Delta t$ is the time step for integration.

An integration scheme for the $NVT$ ensemble using the Nos\'{e}-Hoover chain can also be formulated using the approach of the time-evolution operator . The total Liouville operator for the equations of motion in the Nos\'{e}-Hoover chain method is \begin{equation} iL = iL_1 + iL_2 + iL_{g1} + iL_{g2} + iL_{Tp} + iL_{Tb}, \end{equation} where $iL_{Tp}$ and $iL_{Tp}$ are the Liouville operators for the Nos\'{e}-Hoover chain coupling to the particles and the box, respectively. They are similar to the operator $iL_{T}$ presented in the last subsection. The other four operators are:

\begin{equation} iL_1 = \sum_{i=1}^N \left[\frac{\vec{p}_i}{m_i} + \frac{\mathbf{p}_g}{W_g} \cdot \vec{r}_i\right] \cdot \frac{\partial }{\partial \vec{r}_i}; \end{equation}

\begin{equation} iL_2 = \sum_{i=1}^N \left[\vec{F}_i - \left(\frac{\mathbf{p}_g}{W_g} + \frac{1}{N_f}\frac{\textmd{Tr}[\mathbf{p}_g]}{W_g} \mathbf{I}\right) \cdot \vec{p}_i \right]\cdot \frac{\partial }{\partial \vec{p}_i}; \end{equation}


\begin{equation} iL_{g1} = \frac{\mathbf{p}_g \cdot \mathbf{h}}{W_g} \cdot \frac{\partial }{\partial \mathbf{h}}; \end{equation}


\begin{equation} iL_{g2} = \left(\det[\mathbf{h}] (\boldsymbol{\sigma}-\boldsymbol{\sigma}^{\rm ext}) + \frac{1}{N_f} \sum_{i=1}^N \frac{\vec{p}_i^2}{m_i} \mathbf{I}\right) \cdot \frac{\partial }{\partial \mathbf{p}_g}. \end{equation}


The total time-evolution operator $e^{iL\Delta t} $ for one step can be factorized using the Trotter theorem (the factorization is not unique): \begin{equation} e^{iL} \approx e^{iL_{Tb}\Delta t/2} e^{iL_{Tp}\Delta t/2} e^{iL_{g2}\Delta t/2} e^{iL_2\Delta t/2} e^{iL_{g1}\Delta t} e^{iL_1\Delta t} e^{iL_2\Delta t/2} e^{iL_{g2}\Delta t/2} e^{iL_{Tp}\Delta t/2} e^{iL_{Tb}\Delta t/2}. \end{equation} The two operators related to the box, $e^{iL_{g1}\Delta t}$ and $e^{iL_{g2}\Delta t/2}$, are simple translations. The operators related to the Nos\'{e}-Hoover chain have been discussed in detail in the previous subsection. The remaining operators related to the particles are more complicated. The operator $e^{iL_{1}\Delta t}$ gives \begin{equation} \vec{r}_i(\Delta t)= \mathbf{O}^T \cdot \mathbf{D} \cdot \mathbf{O} \cdot \vec{r}_i(0) + \Delta t \mathbf{O}^T \cdot \tilde{\mathbf{D}} \cdot \mathbf{O} \cdot \vec{v}_i; \end{equation} The operator $e^{iL_{2}\Delta t/2}$ gives \begin{equation} \vec{v}_i(\Delta t/2)= \mathbf{O}^T \cdot \mathbf{\Delta} \cdot \mathbf{O} \cdot \vec{v}_i(0) + \frac{\Delta t}{2m_i} \mathbf{O}^T \cdot \tilde{\mathbf{\Delta}} \cdot \mathbf{O} \cdot \vec{F}_i. \end{equation} Here, the operator $\mathbf{O}$ diagonalizes the box velocity: \begin{equation} \left( \begin{array}{ccc} \lambda_1 & 0 & 0\\ 0 & \lambda_2 & 0\\ 0 & 0 & \lambda_3\\ \end{array} \right) = \mathbf{O} \cdot \frac{\mathbf{p}_g}{W_g} \cdot \mathbf{O}^T, \end{equation} where $\lambda_1$, $\lambda_2$, and $\lambda_3$ are the eigenvalues of the velocity matrix. The other matrices in the above equations are all diagonal with the elements: \begin{equation} D_{\alpha\alpha} = e^{\lambda_{\alpha} \Delta t}; \end{equation} \begin{equation} \tilde{D}_{\alpha\alpha} = e^{\lambda_{\alpha} \Delta t/2} \frac{\sinh\left[\lambda_{\alpha} \Delta t/2\right]}{\lambda_{\alpha} \Delta t/2}; \end{equation}

\begin{equation} \Delta_{\alpha\alpha} = \exp\left[-\left(\lambda_{\alpha} + \frac{1}{N_f} \frac{\textmd{Tr}[\mathbf{p}_g]}{W_g}\right)\frac{\Delta t}{2}\right] ; \end{equation}

\begin{equation} \tilde{\Delta}_{\alpha\alpha} = \exp\left[-\left(\lambda_{\alpha} + \frac{1}{N_f} \frac{\textmd{Tr}[\mathbf{p}_g]}{W_g}\right)\frac{\Delta t}{4}\right] \frac{\sinh\left[\left(\lambda_{\alpha} + \frac{1}{N_f} \frac{\textmd{Tr}[\mathbf{p}_g]}{W_g}\right)\frac{\Delta t}{4}\right]}{\left(\lambda_{\alpha} + \frac{1}{N_f} \frac{\textmd{Tr}[\mathbf{p}_g]}{W_g}\right)\frac{\Delta t}{4}}. \end{equation} As the argument in the sinh function here is small, it is better to evaluate the ratio $\sinh(x)/x$ using a power series expansion: \begin{equation} \frac{\sinh(x)}{x} \approx \sum_{n=0}^m\frac{1}{(2n+1)!} x^{2n}. \end{equation} where $m=5$ is accurate enough.

相互作用势模型(未修订,请不要看)

经验势函数的一般性质

一个多粒子系统的势能可以写为 \begin{equation} U=U\left(\{\vec{r}_i\}_{i=1}^N\right). \end{equation} 该势能应该在欧几里得群的操作下保持不变。据此可以证明系统势能仅仅依赖于两两粒子对之间的距离: \begin{equation} U=U\left((\{r_{ij}\}_{i<j}\right). \end{equation} 如果可以进一步将体系的势能写成 \begin{equation} \label{equation:U-two-body} U= \frac{1}{2}\sum_{i}\sum_{j \neq i} U_{ij} \left(r_{ij} \right), \end{equation} 其中, \begin{equation} U_{ij} \left(r_{ij} \right) = U_{ji} \left(r_{ji} \right), \end{equation} 那么我们称该系统的相互作用势能为两体势(two-body potential)。其中,[math]U_{ij} \left(r_{ij} \right)[/math] 代表粒子 [math]i[/math][math]j[/math] 之间的相互作用势能,仅仅依赖于两粒子的相对距离 [math]r_{ij}[/math]。 两体势系统的势能也可以写成如下等价的形式: \begin{equation} U= \sum_{i}\sum_{j > i} U_{ij} \left(r_{ij} \right). \end{equation} 如果一个体系的势能无法写成以上形式,那么我们称该势能为多体势(many-body potential)。

两体势

Lennard-Jones 势

我们已经在上一章讨论过LJ势。

\begin{equation} U_{ij} = 4 \epsilon \left( \frac{ \sigma^{12} }{ r_{ij}^{12} } - \frac{\sigma^{6} }{ r_{ij}^{6} } \right). \end{equation}

Buckingham-库伦势

该势函数也叫做刚性离子(rigid-ion)势,它由 Buckingham 势 \begin{equation} U_{ij} = A_{ij} \exp\left( - b_{ij} r_{ij} \right) - \frac{C_{ij} }{ r_{ij}^{6} } \end{equation} 和库伦势 \begin{equation} U_{ij} = \frac{q_iq_j}{4\pi\epsilon_0r_{ij}} \end{equation} 组成。我们将在后面的一节专门讨论库伦势。


库伦势

Ewald求和

PPPM等方法

Wolf等方法

多体势

多体势中力的表达式

  • 作者在 2015 年推导出如下力的表达式:

$$ \vec{F}_{i} = \sum_{j \neq i} \vec{F}_{ij}; $$ $$ \vec{F}_{ij} = - \vec{F}_{ji} = \frac{\partial U_{i}}{\partial \vec{r}_{ij}} - \frac{\partial U_{j}}{\partial \vec{r}_{ji}} = \frac{\partial \left(U_{i} + U_{j}\right) }{\partial \vec{r}_{ij}}. $$ 这里, $$ \partial U_{i}/\partial \vec{r}_{ij} = \partial U_{i}/\partial x_{ij} \vec{e}_x + \partial U_{i}/\partial y_{ij} \vec{e}_y + \partial U_{i}/\partial z_{ij} \vec{e}_z $$


集团势(cluster potential)

SW势

键序势(bond-order potential)

键序势也叫 Abell-Tersoff-Brenner势,是 Abell、Tersoff 和 Brenner 等人提出并发展的一类势函数。


Tersoff 势

Tersoff 势有几个稍有不同的变体。我这里介绍 Tersoff 在 1989 年发表的一篇文章中使用的形式,简称为 Tersoff-1989 势。为简单起见,我们考虑单种元素的势函数。

粒子 [math]i[/math] 的势能可以写为: \begin{equation} U_i = \frac{1}{2} \sum_{j \neq i} f_C(r_{ij}) \left[ f_R(r_{ij}) - b_{ij} f_A(r_{ij}) \right]. \end{equation} 其中,[math]f_{C}[/math] 是一个截断函数,当 [math]r_{ij}\lt R[/math] 时取值为 1,当 [math]r_{ij}\gt S[/math] 时取值为 0。在这之间,该函数为 \begin{equation} f_{C}(r_{ij}) = \frac{1}{2} \left[ 1 + \cos \left( \pi \frac{r_{ij} - R_{ij}}{S_{ij} - R_{ij}} \right) \right]. \end{equation}

排斥函数 [math]f_{R}[/math] 和吸引函数 [math]f_{A}[/math] 为 \begin{equation} f_{R}(r) = A e^{-\lambda r_{ij}}; \end{equation} \begin{equation} f_{A}(r) = B e^{-\mu r_{ij}}. \end{equation}

键序为 \begin{equation} \label{equation:bij} b_{ij} = \left(1 + \beta^{n} \zeta^{n}_{ij}\right)^{-\frac{1}{2n}}, \end{equation} 其中, \begin{equation} \zeta_{ij} = \sum_{k\neq i, j}f_C(r_{ik}) g_{ijk}, \end{equation} \begin{equation} g_{ijk} = 1 + \frac{c^2}{d^2} - \frac{c^2}{d^2+(h-\cos\theta_{ijk})^2}. \end{equation}

在以上表达式中,有如下参数:[math]A[/math], [math]B[/math], [math]\lambda[/math], [math]\mu[/math], [math]\beta[/math], [math]n[/math], [math]c[/math], [math]d[/math], [math]h[/math], [math]R[/math], [math]S[/math]

两体泛函势

EAM 势

  • 原子 [math]i[/math] 的势能为

\begin{equation} U_i = \frac{1}{2} \sum_{j\neq i} \phi(r_{ij}) + F (\rho_i). \end{equation} 这里,含有 [math]\phi(r_{ij})[/math] 的部分是两体势,[math]F(\rho_i)[/math] 即为嵌入势。嵌入势是 [math]i[/math] 粒子处电子密度 [math]\rho_i[/math] 的函数。粒子 [math]i[/math] 所在点的电子密度是由它的邻居贡献的: \begin{equation} \rho_i = \sum_{j\neq i} f(r_{ij}). \end{equation}


  • 练习。推导如下表达式:

\begin{equation} \frac{\partial U_i}{\partial \vec{r}_{ij}} = \frac{1}{2} \phi'(r_{ij}) \frac{\partial r_{ij}} {\partial \vec{r}_{ij}} + F'(\rho_i) f'(r_{ij}) \frac{\partial r_{ij}} {\partial \vec{r}_{ij}}. \end{equation}


推广的 EAM 势

介绍 MEAM 势。


拓扑力场

OPLS 等

反应力场

reaxFF 等

机器学习势

平衡态分子动力学模拟(未修订,请不要看)

结构性质

状态方程

径向分布函数

弹性常数

热力学性质

简单的热力学性质

热力学响应函数

振动态密度与相关热力学性质

基于线性响应理论的模拟方法

时间关联函数

  • 在统计力学中,设有两个依赖于时间的物理量 [math]A(t)[/math][math]B(t)[/math],我们定义这两个量之间的时间关联函数 (time correlation function) [math]C(t)[/math] 为:

\begin{equation} C(t) = \left\langle A(t_0) B(t_0+t) \right\rangle. \end{equation} 对这个公式的说明如下:

    • 如果两个物理量相同,[math]A=B[/math],那么上式代表物理量 [math]A(t)[/math]自关联函数 (auto-correlation function)。
    • 上式中~$t$ 是某个时间间隔,叫做关联时间 (correlation time),而时间关联函数~$C(t)$ 是关联时间~$t$ 的函数。
    • 尖括号在统计物理中代表系综平均,但在分子动力学模拟中一般代表“时间”平均,其中“时间”指的不是上述“关联时间”~$t$,而是“时间原点”~$t_0$。一般将~$t_0$ 写作~0,即常写~$C(t) = \left\langle A(0) B(t) \right\rangle$。
    • 我们这里考虑的是平衡系统,即在每一个时间原点~$t_0$,系统都处于平衡态。所以,不同的时间原点在物理上是等价的,从而可以对时间原点求平均。


在分子动力学模拟中,我们只能得到一条离散的相轨迹(phase trajectory),故在实际计算中,与尖括号对应的时间平均将由求和方式表示。

首先,我们假设先在控制温度的情况下让系统达到平衡。然后在不控制温度的情况下(即在微正则系综)模拟了~$N_p$ 步(下标~$p$ 是~production 的意思),步长为~$\Delta t$。与计算静态热力学量的情形类似,我们不需要将每一步的数据都保存(因为相邻步的数据有关联性,保存得过频无益)。我们假设每~$N_s$ (下标~$s$ 是~sampling 的意思) 步保存一次数据,并称~$N_s$ 为取样间隔~(sampling interval)。我们假设~$N_p$ 是~$N_s$ 的整数倍,记 \begin{equation} \frac{N_p}{N_s} = N_d \end{equation} 代表记录数据的总步数(下标~$d$ 是~data 的意思),并记 \begin{equation} N_s \Delta t = \Delta \tau. \end{equation}

根据上面的讨论与记号,我们可以将关联时间为~$t$ 的关联函数表达为: \begin{equation} C(t) = \frac{1}{M}\sum_{m=1}^{M} A(m \Delta \tau) B(m\Delta \tau+t). \end{equation} 其中,$M$ 是求平均时用的时间原点数目。原则上,对不同的关联时间~$t$,可被利用的时间原点的最大数目~$M$ 是不同的(注意:$t$ 一定是~$\Delta \tau$ 的倍数): \begin{equation} M = N_d - \frac{t}{\Delta \tau}. \end{equation} 关联时间等于零时,最大时间原点数目~$M=N_d$;关联时间等于几个时间步长,$M$ 的值就要减少几。在实际的模拟中,要考虑的最大关联数据量~$N_c$ (下标~$c$ 是~correlation 的意思) 通常远小于~$N_d$,如果希望得到比较精确的模拟结果的话。所以,为了编程方便(特别是用并行计算时),通常将~$M$ 的值统一取为 \begin{equation} M = N_d - N_c. \end{equation} 这样的话,虽然浪费了一小部分数据,但好处是与每个关联时间对应的关联函数值的统计平均采用了统一的数据量(即时间原点数)。从统计的角度来看,这样做是更合理的。

我举一个非常简单的例子。假设保存了~$N_d=100$ 组数据,并且要计算~$N_c=10$ 个关联函数值。首先确定~$M=100-10=90$。这~10 个关联函数值为 \begin{equation} C(0) = \frac{1}{90}\sum_{m=1}^{90} A(m \Delta \tau) B(m\Delta \tau); \end{equation} \begin{equation} C(\Delta \tau) = \frac{1}{90}\sum_{m=1}^{90} A(m \Delta \tau) B(m\Delta \tau + \Delta \tau); \end{equation} \begin{equation} C(2\Delta \tau) = \frac{1}{90}\sum_{m=1}^{90} A(m \Delta \tau) B(m\Delta \tau + 2\Delta \tau); \end{equation} \begin{equation} C(9\Delta \tau) = \frac{1}{90}\sum_{m=1}^{90} A(m \Delta \tau) B(m\Delta \tau + 9\Delta \tau). \end{equation}

Green-Kubo 理论

格林-久保公式实际上是一类公式,它们将非平衡过程的输运系数与平衡态中相应物理量的涨落相联系。格林-久保公式是说,输运系数等于自关联函数对关联时间的积分。例如,扩散系数是速度自关联函数的积分;粘性系数是压力自关联函数的积分;热导率是热流自关联函数的积分、等等。


速度自关联、均方位移与扩散系数

Velocity autocorrelation (VAC) is an important quantity in MD simulations. On the one hand, its integral with respect to the correlation time gives the running diffusion constant, which is equivalent to that obtained by a time derivative of the mean square displacement (MSD). On the other hand, its Fourier transform is the phonon density of states (PDOS).

\subsection{Running diffusion coefficient}

The VAC is a single-particle correlation function. This means that we can define the VAC for individual particles. For particle $i$, the VAC along the $x$ direction is defined as \begin{equation} \langle v_{xi}(0) v_{xi}(t) \rangle. \end{equation} Then, one can define the mean VAC for any number of particles. In the current version of GPUMD, it is assumed that one wants to calculate the mean VAC in the whole simulated system: \begin{equation} \boxed{ \text{VAC}_{xx}(t) = \frac{1}{N} \sum_{i=1}^{N} \langle v_{xi}(0) v_{xi}(t) \rangle }. \end{equation} The order between the time-average (denoted by $\langle \rangle$) and the space-average (the average over the particles) can be changed: \begin{equation} \boxed{ \text{VAC}_{xx}(t) = \left\langle \frac{1}{N} \sum_{i=1}^{N} v_{xi}(0) v_{xi}(t) \right\rangle }. \end{equation} Using the same conventions as in the case of HAC calculations, we have the following explicit expression for the VAC: \begin{equation} \label{equation:VAC} \text{VAC}_{xx}(n_c\Delta \tau) = \frac{1}{(N_d-N_c)N} \sum_{m=0}^{N_d-N_c-1} \sum_{i=1}^{N} v_{xi}(m\Delta \tau) v_{xi}((m+n_c)\Delta \tau), \end{equation} where $n_c = 0, 1, 2, \cdots, N_c-1$. The algorithm for calculating the VAC is quite similar to that for calculating the HAC and it thus omitted.

After obtaining the VAC, we can calculate the running diffusion constant $D_{xx}(t)$ as \begin{equation} \boxed{ D_{xx}(t) = \int_0^{t} dt' ~\text{VAC}_{xx}(t') }. \end{equation} One can prove that this is equivalent to the time-derivative of the MSD, i.e., the Einstein formula: \begin{equation} \boxed{ D_{xx}(t) = \frac{1}{2} \frac{d}{dt} \Delta x^2(t) }, \end{equation} where the MSD $\Delta x^2(t)$ is defined as \begin{equation} \boxed{ \Delta x^2(t) = \left\langle \frac{1}{N} \sum_{i=1}^{N} \left[ x_i(t) - x_i(0) \right]^2 \right\rangle = \frac{1}{N} \sum_{i=1}^{N} \left\langle

 \left[ x_i(t) - x_i(0) \right]^2

\right\rangle }. \end{equation}

Here is the proof. Starting from the relation between position and velocity, \begin{equation} x_i(t) - x_i(0) = \int_{0}^{t}dt' v_{xi}(t'), \end{equation} we have \begin{equation} [x_i(t) - x_i(0)]^2 = \int_{0}^{t} dt' v_{xi}(t') \int_{0}^{t}dt v_{xi}(t)= \int_{0}^{t} dt' \int_{0}^{t}dt v_{xi}(t') v_{xi}(t). \end{equation} Then, the MSD can be expressed as \begin{equation} \Delta x^2(t) = \frac{1}{N}\sum_{i=1}^{N} \int_{0}^{t} dt' \int_{0}^{t}dt \left\langle v_{xi}(t') v_{xi}(t) \right\rangle. \end{equation} Using Lebniz's rule, we have \begin{equation} D_{xx}(t) = \frac{1}{2} \frac{d}{dt} \Delta x^2(t) = \frac{1}{N}\sum_{i=1}^{N} \int_{0}^{t} dt' \left\langle v_{xi}(t) v_{xi}(t') \right\rangle, \end{equation} which can be rewritten as \begin{equation} D_{xx}(t) = \frac{1}{N}\sum_{i=1}^{N} \int_{0}^{t} dt' \left\langle v_{xi}(0) v_{xi}(t'-t) \right\rangle. \end{equation} Letting $\tau=t'-t$, we get (note that here $t$ is considered as a constant) \begin{equation} D_{xx}(t) = \frac{1}{N}\sum_{i=1}^{N} \int_{-t}^{0} d\tau \left\langle v_{xi}(0) v_{xi}(\tau) \right\rangle, \end{equation} which can be rewritten as \begin{equation} D_{xx}(t) = \frac{1}{N}\sum_{i=1}^{N} \int_{-t}^{0} d\tau \left\langle v_{xi}(-\tau) v_{xi}(0) \right\rangle. \end{equation} Letting $t'=-\tau$, we finally get \begin{equation} D_{xx}(t) = \frac{1}{N}\sum_{i=1}^{N} \int_{0}^{t} dt' \left\langle v_{xi}(t') v_{xi}(0) \right\rangle =\int_0^t dt' ~\text{VAC}_{xx}(t'). \end{equation} We thus have derived the Green-Kubo formula from the Einstein formula.

In summary, \begin{itemize} \item The derivative of half of the MSD gives the running diffusion coefficient. \item The integral of the VAC gives the running diffusion coefficient. \item One can obtain the MSD by integrating the VAC twice (numerically). \end{itemize}


It is interesting that the same VAC can be used to compute the PDOS, as first demonstrated by Dickey and Paskin. The PDOS is simply the Fourier transform of the normalized VAC: \begin{equation} \rho_x(\omega) = \int_{-\infty}^{\infty} dt e^{i\omega t}~\text{VAC}_{xx}(t). \end{equation} Here, $\text{VAC}_{xx}(t)$ should be understood as the normalized function $\text{VAC}_{xx}(t)/\text{VAC}_{xx}(0)$. Although it looks simple, it does not mean that one can get the correct PDOS by a naive fast Fourier transform (FFT) routine. Actually, this computation is very cheap and we do not need FFT at all. What we need is a discrete cosine transform. To see this, we first note that, by definition, $\text{VAC}_{xx}(-t) = \text{VAC}_{xx}(t)$. Using this, we have \begin{equation} \rho_x(\omega) = \int_{-\infty}^{\infty} dt \cos (\omega t)~\text{VAC}_{xx}(t). \end{equation} Because we only have the VAC data at the $N_c$ discrete time points, the above integral is approximated by the following discrete cosine transform: \begin{equation} \rho_x(\omega) \approx \sum_{n_c=0}^{N_c-1} (2-\delta_{n_c0}) \Delta \tau \cos (\omega n_c \Delta \tau)~\text{VAC}_{xx}(n_c \Delta \tau). \end{equation} Here, $\delta_{n_c0}$ is the Kronecker $\delta$ function and the factor $(2-\delta_{n_c0})$ accounts for the fact that there is only one point for $t = 0$ and there are two equivalent points for $t \neq 0$. Last, we note that a window function is needed to suppress the unwanted Gibbs oscillation in the calculated PDOS. In GPUMD, the Hann window $H(n_c)$ is applied: \begin{equation} \rho_x(\omega) \approx \sum_{n_c=0}^{N_c-1} (2-\delta_{n_c0}) \Delta \tau \cos (\omega n_c \Delta \tau)~\text{VAC}_{xx}(n_c \Delta \tau) H(n_c); \end{equation} \begin{equation} H(n_c) = \frac{1}{2} \left[ \cos \left( \frac{\pi n_c}{N_c} \right) + 1 \right]. \end{equation}

Here are some comments on the normalization of the PDOS. In the literature, one usually uses an arbitrary unit for the PDOS, but it actually has a dimension of [time], and an appropriate unit for it can be 1/THz or ps. The normalization of $\rho_x(\omega)$ can be determined by the inverse Fourier transform: \begin{equation} \text{VAC}_{xx}(t)

= \int_{-\infty}^{\infty} \frac{d\omega}{2\pi} e^{-i\omega t}\rho_x(\omega).

\end{equation} As we have normalized the VAC, we have \begin{equation} 1 = \text{VAC}_{xx}(0)

= \int_{-\infty}^{\infty}
\frac{d\omega}{2\pi}\rho_x(\omega).

\end{equation} Because $\rho_x(-\omega)=\rho_x(\omega)$, we have \begin{equation} \int_{0}^{\infty}

\frac{d\omega}{2\pi}\rho_x(\omega) = \frac{1}{2}.

\end{equation} The calculated PDOS should meet this normalization condition (approximately).


热流自关联与热导率

根据热力学第二定律,一个系统在不受外场作用时,若其内部有热力学性质的不均匀性,则它一定处于非平衡的状态,并有向平衡态靠近的趋势。这种由热力学性质的不均匀性导致的热力学过程叫做输运过程~(transport process),相应的现象叫做输运现象。例如,温度的不均匀性导致能量的输运(热传导现象);粒子数密度的不均匀性导致粒子的输运(扩散现象)。将一个系统置于两个温度不同的热源之间,最终会在系统内建立一个稳定的~(不随时间变化的) 温度分布。我们说这样的系统处于一个稳态~(steady state),但不处于一个平衡态~(equilibrium state)。稳态和平衡态都是不依赖与时间的,但前者属于非平衡态。

上述不均匀性都是由相应的不均匀的物理量的梯度来量化的。这里,我们只研究热输运,而且假设输运方向沿着一个特定方向(假设是~$x$ 方向)的情形。热传导现象的宏观规律由傅里叶定律描述。 傅里叶定律是说热流密度~(heat flux,或者~heat current) $J$,即单位时间穿过单位面积的热量,在数量上正比于温度梯度~$\frac{dT}{dx}$: \begin{equation} J = - \kappa \frac{dT}{dx}. \end{equation} 这里的~$\kappa$ 就反映了热量输运的难易程度:$\kappa$ 越大代表热量越容易被输运。这样的物理量被称为输运系数~( transport coefficient)。具体到热传导,输运系数~$\kappa$ 叫做热导率~(thermal conductivity)。注意等式右边有个负号,它表示热量的传导方向与温度梯度的方向相反,指向温度降低的方向~(一个物理量的梯度的方向指向它增加的方向)。在国际单位制中,温度梯度的单位为~K/m,热流密度的单位是~W/m$^2$,故热导率的单位是~$\text{W}\text{m}^{-1}\text{K}^{-1}$。


对热导率的计算有如下格林-久保公式: \begin{equation} \kappa_{\mu\nu}(t) = \frac{V}{k_B T^2} \int_0^{t} dt' C_{\mu\nu} (t'). \end{equation} 其中,$\kappa_{\mu\nu}(t)$ ($\mu, \nu = x, y, z$) 是热导率张量,$t'$ 是关联时间,$k_B$ 是~Boltzmann 常数, $T$ 是温度, $V$ 是体积,$C_{\mu\nu}(t)$ 是热流自关联函数(heat current autocorrelation function,常简称为~HCACF)。上式计算的跑动热导率(running thermal conductivity)。关于跑动热导率的技术细节,将在具体的范例中讨论。


热流自关联函数的表达式如下: \begin{equation}

C_{\mu\nu}(t) = \langle J_{\mu}(0) J_{\nu}(t) \rangle.

\end{equation} 其中,尖括号表示统计平均,在分子动力学模拟中指对时间原点的平均,$J_{\mu}$ ($\mu = x, y, z$) 是热流。下一节讨论关联函数;下下节讨论热流的具体表达式。


如果研究的是三维的各向同性(isotropic)的系统,则热导率张量的非对角分量一定为零,并可将最终计算的热导率取为对角分量的平均值: \begin{equation} \kappa = \frac{\kappa_{xx} + \kappa_{yy} + \kappa_{zz}}{3}. \end{equation} 使用格林-久保方法时要注意边界条件的选取: \begin{itemize} \item 如果模拟的是三维块体(bulk)系统,则每个方向都要使用周期边界条件。 \item 如果要研究准两维系统(如薄膜或者两维材料),则在垂直于薄膜的方向用自由边界条件,在平行于薄膜的方向用周期边界条件。而且,此时垂直方向热导率的计算结果无意义。 \item 如果要研究的是准一维系统(如纳米线或者纳米管),则在垂直于线或者管的方向都要用自由边界条件,在平行于线或者管的方向用周期边界条件。而且,此时只有平行方向的热导率结果才有意义。 \end{itemize}


热流定义为能量密度矩的时间导数: \begin{equation} \vec{J} = \frac{1}{V} \frac{d}{d t} \sum_i \vec{r}_i E_i. \end{equation} 其中, \begin{equation} E_i = \frac{1}{2} m_i \vec{v}_i^2 + U_i \end{equation} 是第~$i$ 个粒子的总能量,$U_i$ 是势能部分,$\frac{1}{2} m_i \vec{v}_i^2$ 是动能部分。

首先,用求导的莱布尼茨法则可得 \begin{equation}

 \vec{J} = \sum_i \vec{v}_i E_i + \sum_i \vec{r}_i \frac{d}{d t} E_i

\end{equation} 通常将上式右边的两项分别称为动能项 \begin{equation} \vec{J}_{\textmd{kin}} = \sum_i \vec{v}_i E_i \end{equation} 和势能项 \begin{equation} \vec{J}_{\textmd{pot}} = \sum_i \vec{r}_i \frac{d}{d t} E_i. \end{equation} 合起来,有 \begin{equation} \vec{J} = \vec{J}_{\textmd{kin}} + \vec{J}_{\textmd{pot}}. \end{equation} 动能项不需要再推导了。利用动能定理 \begin{equation} \frac{d}{dt}\left(\frac{1}{2}m_i\vec{v}_i^2\right)= \vec{F}_i \cdot \vec{v}_i, \end{equation} 可以将势能项写为 \begin{equation}

 \vec{J}_{\textmd{pot}} = \sum_i \vec{r}_i (\vec{F}_i \cdot \vec{v}_i)
 + \sum_i \vec{r}_i \frac{d U_{i}}{d t}.

\end{equation}


在两体势(two-body potentials)情形中,系统总势能可表达为 \begin{equation}

 U =  \frac{1}{2} \sum_{i}  \sum_{j \neq i} U_{ij}(r_{ij}).

\end{equation} 可以推导如下力的表达式: \begin{equation}

 \vec{F}_i = \sum_{j \neq i} \vec{F}_{ij},

\end{equation} \begin{equation}

 \vec{F}_{ij}
 = \frac{\partial U_{ij}}{\partial \vec{r}_{ij}}
 = - \vec{F}_{ji}.

\end{equation} 其中,$\vec{F}_{ij}$ 是第~$i$ 个粒子受到的来自于第~$j$ 个粒子的力, \begin{equation}

\vec{r}_{ij} \equiv \vec{r}_{j} - \vec{r}_{i},

\end{equation} 是从第~$i$ 个粒子指向第~$j$ 个粒子的位置差矢量。

\vspace{0.5cm} \hrule \hrule 练习。定义 \begin{equation} U_i = \frac{1}{2}\sum_{j\neq i} U_{ij}, \end{equation} 推导如下公式 \begin{equation}

 \vec{J}_{\textmd{pot}}
 = \frac{1}{2}\sum_i \sum_{j \neq i}
   \vec{r}_{i} [\vec{F}_{ij} \cdot (\vec{v}_i + \vec{v}_j)].

\end{equation} 在分子动力学模拟中,如果采用了周期边界条件,一个宏观物理量就不可能依赖于绝对坐标~$\vec{r}_i$, 而应该依赖于相对坐标。证明上面的热流表达式等价于 \begin{equation}

 \vec{J}_{\textmd{pot}}
 = - \frac{1}{4}\sum_i \sum_{j \neq i}
   \vec{r}_{ij} [\vec{F}_{ij} \cdot (\vec{v}_i + \vec{v}_j)].

\end{equation} 再证明上式等价于下面的不怎么对称的形式: \begin{equation}

 \vec{J}_{\textmd{pot}}
 = - \frac{1}{2}\sum_i \sum_{j \neq i}
   \vec{r}_{ij} [\vec{F}_{ij} \cdot \vec{v}_i].

\end{equation} \hrule \hrule \vspace{0.5cm}


在分子动力学模拟中,位力~(virial) 张量定义为 \begin{equation}

\textbf{W} = \sum_i \textbf{W}_{i},

\end{equation} 其中,$\textbf{W}_i$ 是单粒子位力: \begin{equation}

\label{equation:per_atom_virial}
\textbf{W}_i
= -\frac{1}{2} \sum_{j \neq i}\vec{r}_{ij} \otimes \vec{F}_{ij}.

\end{equation} 于是有: \begin{equation} \label{equation:j_pot_pair_stress}

 \vec{J}_{\textmd{pot}} = \sum_{i} \textbf{W}_{i} \cdot \vec{v}_i.

\end{equation} 这个公式就是~LAMMPS 中用的热流公式,但它只适用于两体势,不适用于多体势(many-body potentials)。对多体势的讨论,请参看我的论文。


压力自关联与粘滞系数

非平衡态分子动力学模拟(未修订,请不要看)

非平衡稳态模拟

热传导模拟

  • 前已提及,热传导是温度梯度引起能量传递的现象,属于非平衡过程。另一种计算热导率的方法就是基于非平衡分子动力学模拟的。在该方法中,首先通过某种方式建立一个温度梯度,并等待足够长的时间,使得系统达到稳态。在达到稳态后,对系统的温度梯度和非平衡稳态热流进行测量,然后就可以根据傅里叶定律来计算热导率。
  • 必须区分热导率 (thermal conductivity) 和热导 (thermal conductance) 这两个不同的概念。考虑横截面积为 [math]A[/math],长度为 [math]d[/math] 的一个长方体材料,热导 [math]K[/math] 与热导率 [math]\kappa[/math] 有如下关系:

\begin{equation} K = \kappa \frac{A}{d}. \end{equation} 更常用的是单位面积的热导,我们用符号 [math]G[/math] 表示。对于上述长方体材料,单位面积的热导为 \begin{equation} G = \frac{K}{A} = \frac{\kappa}{d}. \end{equation} 这样定义的 [math]G[/math] 的单位为 Wm-2K-1。为简单起见,我们以后都简称 [math]G[/math] 为热导。我们也将简称热流密度为热流。这些都是文献中常用的简称,虽然不严谨。

  • 产生温度梯度的方法有多种。在每一种方法中,都将模拟系统在热传导方向分割为若干块

(blocks)。每个块的大小可以不等,但一般为了处理数据的方便,取为相等。将其中的两块分别设置为热源(heat source)区域和热汇(heat sink)区域。热源和热汇区域的位置选取与模拟的边界条件相关,具体如下:

    • 如果采用周期边界条件,则块的个数一般取为偶数~$2N$(假设从~1 到~$2N$ 进行标记),并可以取第~$1$ 块为热源,取第~$N+1$ 块为热汇,或者反过来。
    • 如果采用固定边界条件,则一般将传导方向的前后两端一小层原子固定,并将固定层附近的块取为热源和热汇。此时块的个数可以不为偶数。


  • 产生温度梯度的方法:
    • 动量交换法。该方法由~M\"{u}ller-Plathe 提出。
    • 速度重标法。这个方法首先被 Ikeshoji 和 Hafskjold 提出,后来又被 Jund 和 Jullien 提出(估计后者不知道前者的工作)。
    • 局部热浴法。这个方法好像不归功于任何个人。我不知道是谁最先用这个方法的,但我认为这个方法是最方便的。在动量交换法中,动量交换的频率这个参数一般要用试错法来确定。交换过频,引起的温度梯度会过大。同样,在速度重标法中,每次给热区增加的能量值也一般要用试错法来确定。相比之下,在局部热浴法中,我们可以直接设置热源和热汇的温度值。例如,如果要模拟一个系统在~300 K 时的热传导,我们可以将热源的温度设置为~310 K,将热汇的温度设置为~290 K,从而直接保证温度梯度具有合理的值。
  • 局部热浴法的算法为:
    • 在系统达到平衡态之后,撤掉整体热浴,分别对热源和热汇区域施加局部热浴,并取合适的温度值。常用的局部热浴有~Nose-Hoover (chain) 热浴和郎之万(Langevin)热浴。
    • 热流大小可以通过系统和热浴交换的能量来计算(能量守恒)。
    • 在系统达到稳态之后,开始测量温度梯度,然后根据傅里叶定律计算热导率。

例子:Kapitza 热阻的计算

  • 远离平衡态模拟
   * 导热
   * 拉伸
   * 摩擦
   * 切削

分子动力学模拟的并行计算(未修订,请不要看)

openMP并行计算

MPI并行计算

CUDA并行计算