在分子模拟中，一个重要的技术是能够维持体系的温度在一定范围内不变（或者严格恒温），来模拟真实的实验条件。这便是我们常说的恒温分子动力学模拟，或在NVT系综下的模拟。在这一模拟中，我们“近似的看作”体系在恒定的温度下演化。

等等，恒温模拟的温度不应该是 恒定（constant）的吗？为什么说是“近似”？恒温模拟到底是怎么实现的？

为了解决这个问题，我们需要回答下面一些问题：

1. 什么是温度？
2. 理论上如何定义“恒温”？实际上如何定义“恒温”？
3. 如何实现实际意义上的“恒温”？
4. 如何评价不同算法的优劣？

要回答着一系列问题，我们需要从一些热力学和统计力学说起。

什么是温度？什么是恒温？

温度就是温度计的示数吗？

有的同学会说，温度（temperature） 这个东西太熟悉了，不就是通过温度计测量物体得到的值嘛。这个说法从某种角度是正确的，这是站在**“宏观热力学（Thermodynamics）”**的角度去看待温度。对于宏观物体来说，一个物体的温度就是一个固定的数值，一个准确的热力学态，一个能够在相图上标出来的点。

比如，在一个标准大气压下（但也有别的说法，如IUPAC变更参考态为标准态的另一个数值，即1 bar下的数值），冰的熔点就是273.15K, 三相点就是273.16K，水的沸点就是371.15K，天王老爷来了也不会变。

那么这里的宏观热力学是什么意思呢，这是因为宏观热力学是巨大量粒子组成的宏观物体的状态（这时你基本可以认为微观粒子个数就是无限的）。它最突出的特点就是正常状态下（即没有相变的时候）系统的涨落很小。大白话说就是你测量的物体的温度不会发生很大改变。注意，这里的涨落不是你的温度计的测量误差，而是假设你的温度计无限准的情况下，物体自身内秉性存在的温度变化。

但这里有一个问题，热力学态的温度，到底是怎么定义的？如果你翻变了经典热力学的书，你可能会发现它本身并没有严格的定义。最早的热力学温度，是通过物质参考态来定义的。即规定冰点为0摄氏度，水的沸点为100摄氏度，中间的温度数值做插值；也可以通过汞的体积随温度线性变化得到另一种参考态温度；也可以通过理想气体方程定义（推导）温度。参考态温度之间可以相互转化，但是这种基于物质参考态的定义并没有触及温度的本质。

那粒子体系怎么插温度计？

那么从宏观热力学的角度出发，在微观尺度上的参考物质和参考温度是啥呢？这时你就会发现，你根本无法在你的粒子体系（比如你的模拟盒子）里插入一个温度计，更不用说去“测量”温度了。

这时，教科书会告诉你，温度反应了粒子运动的动能。那是不是只要知道了这个粒子的动能，就知道了这个粒子的温度？答案：不是！！！

非常关键的一件事情：

温度是定义在一个粒子的集合上的。单一粒子不存在“温度”，它的温度没有意义。

举个例子，现在你有一个微观粒子体系，那么温度是在给定的微观态下，体系中所有粒子的动能平均值的表达式。根据能均分定理：

因为分子动力学大部分情况下遵守波恩奥本海默近似，因此下文的所有统计力学表达式均为经典统计力学。

$$k_{B}T=m⟨v_{\alpha }^2⟩=\frac{⟨p_{\alpha }^2⟩}{2m}$$

这里 $k_B$是Boltzmann常数，$T$是微观态的温度，$m$是粒子的质量，$v_{\alpha }$是这个粒子第$\alpha$个自由度分量的速度；也可以等价的用动量$p$表示。因此粒子的动量/速度就是微观体系的表达式。

（理论上）什么是恒温？（Real）

在宏观上很好理解，根据热力学定律，一个宏观物体在热平衡（即与环境温度相等或不存在热传导）下，温度就是一个恒定的数值。自然是“恒温”。

那么由此推出的结论：微观态下粒子的动能保持不变，就是“恒温“。这样对吗？

显然是不对的。因为我们都知道粒子存在热运动，热运动是一种随机的过程，这种随机性必然会改变粒子的运动方式，也就是改变粒子的动量（动量是矢量，既可以改变大小，也可以改变方向）。

考虑到微观下，温度是体系粒子的平均性质，那么我们换一种说法：”微观体系的粒子平均动能保持不变，就是恒温“。这样就更加严谨了。

（模拟中）什么是恒温？（Practical）

等下，我们已经得到了微观体系的温度定义。那么我只要求一下模拟盒子里面的粒子的平均动量，并保持它不变不就好了吗？如果我的盒子有N个粒子，那就有3N个自由度，这个表达式就会是：

$$\hat{T}=\sum _{i=0}^N\sum _{\alpha =1}^3\frac{{p_{\alpha }^i}^2}{2k_B{m}_i}$$

但请仔细想一下，这个表达式和上面的公式真的想等吗？回想我们在概率统计中如何表达一个物理量的期望值，它应当是概率密度与物理量函数的函数内积，也就是说（假设在正则系综下，概率密度满足$P(p,q)\propto e^{-\beta H(p,q)}$）: $$T=\frac{1}{Z}\int ...\int e^{-\frac{H(\text{p},\text{q})}{k_{B}T}}\frac{{\text{p}}^2}{2k_{B}m}\mathrm{d}\text{p}\mathrm{d}\text{q}\approx \sum _{i=0}^N\sum _{\alpha =1}^3\frac{{p_{\alpha }^i}^2}{2k_B{m}_i}=\hat{T}$$

没错，两个表达式并不是完全相等的。因为我们在实际计算的时候，我们不可能模拟无限的个粒子，因此左侧的积分永远不可能真正实现。我们的盒子只有有限个粒子，而且粒子也并不能连续地覆盖整个相空间（即由动量和位置张成的空间）。所以我们只能得到一个有限粒子(finite particles)的非遍历态(non-ergodic)系统。

模拟恒温 和 理论恒温 的区别 —— 盲人摸象

以下的有限系统实验均假设有限时间采样，即不能遍历所有微观态。

！最！大！区！别！：在有限粒子的有限模拟下，温度是浮动的。因为即便系综中粒子动量的均值是恒定的，微观粒子体系也存在局部涨落。这一现象可以通过涨落理论来解释。涨落可以简单理解为这个物理量在给定系综（或给定宏观态）下的微观方差。

温度的涨落

在概率统计中，一个物理量的方差$E[(X-E(X{)}^2)]=E[X^2]-E[X{]}^2$. 在正则系综下（即NVT系综；恒温系综）, $P(\text{p})\propto e^{p^2/2m}$:

$$⟨p^2⟩=\int \mathrm{d}\text{p}{p}^{2}P(p)=\frac{3m}{\beta }$$ $$⟨p^4⟩=\int \mathrm{d}\text{p}{p}^{4}P(p)=15{\left(\frac{m}{\beta }\right)}^2$$

$$\frac{{\sigma }_{p^2}^2}{{⟨p^2⟩}^2}=\frac{2}{3}$$

Since we only have finite particles, let's substitue the temperature of finite system $\hat{T}$ into it:

$$\frac{{\sigma }_T^2}{{⟨T⟩}^2}=\frac{2}{3N}$$

所以，温度的涨落个体系的粒子的个数是有关的。当$N$不是无穷$\inf$的时候，温度就回存在局部的涨落。你可以想象，一个无限体系中的一小部分的温度存在方差，但是当无限个小部分平均在一起，总体均值还是不变的。

这可以类似于盲人摸象。你摸到的那一块（有限粒子体系）只是系综（无限粒子）的一部分，这一小部分可以存在局部的温度涨落。只有把整只大象拼起来，才能得到体系的温度。

因此，在模拟中，你会发现，即便你指定了NVT模拟，体系的温度也是在指定温度的附近变化，而不是恒定在一个值。

Fun Fact: 把温度恒定在一个值的分子模拟算法不能得到正确的正则系综分布。

如何实现（模拟上的）恒温？

虽然我们永远无法得知系综的温度，但也不要灰心。因为在各态历经假设(ergodicity)下，只要我们模拟时间足够长，体系温度的时间平均值与系综的真实温度会不断的减小，直到小于我们预期的误差。

那么，如何在模拟中实现控制体系的温度呢？

Velocity scaling won't work

最简单的想法，那我每一步改变粒子的动量，让它乘以系数后变成我想要的温度不就可以了吗？当然不行，上面我们讲过，温度存在局部的涨落，有限体系的温度内秉性的浮动，使用scaling的方法从理论上就破坏了正则系综分布。事实上，很多数值试验已经表明，这种velocity-scaling的办法会导致错误的系综分布。因为理论上的证明已经很清楚，我们不再进行数值试验。

What's your goal to NVT?

说到这里，其实所谓的“恒温模拟”的主要目标是实现在NVT正则系综下的采样。我们可能并不关心真实的模拟轨迹是什么样子的(区别于传统的Numerical analysis)，而更关心采样分布是否正确，例如是否符合玻尔兹曼分布。

实际上我们大致有两种思路实现NVT系综采样：

• 随机温度控制
• 全局温度控制

随机控温 —— 模拟粒子随机碰撞 (Stochastic Thermostats)

进一步讲，考虑到温度是由粒子的热运动产生的，而热量是由粒子的碰撞传递的。那么如果我们可以通过把体系与一个热浴（heat bath）耦合，实现热量的传递，岂不是就可以稳定系统的温度了。没错！下一步就是如何考虑模拟这种随机碰撞。

• Andersen 热浴 Andersen Thermostat: Andersen 想到，在模拟的每一步后，可以随机地把一些粒子的速度改变，使粒子的速度在时间尺度上符合Maxwell-Boltzmann分布，就可以实现对粒子速度的控制。
• Langevin 热浴 Langevin Dynamics: 从粒子热运动的角度考虑，布朗运动是最能刻画这一行为的方程。而为了考虑到热量的吸收与耗散，我们需要在随机方程中添加摩擦向来代表热量向外传递。这就是Langevin Dynamics的基本思想

全局控温 -- 拓展的拉格朗日力学 (Extended Lagrangian Dynamics)

这一类热浴没有任何随机过程，而是想办法去改进velocity scaling的方式，使得速度能够存在涨落，从而保证采样的正确性。从Hamiltonian的观点出发，在Hamiltonian or Lagrangian Dynamics 下，系统的Hamiltonian是守恒的（即能量守恒）。但NVT系综下，体系和环境存在热交换，体系的能量是不守恒的。如果我们需要Hamiltonian能够产生变化，就需要给体系一些额外的自由度 —— 这就是Extended Lagrangian Dynamics 的基本思想。

• Nose Thermostat and Nose-Hoover Chain. Nose从Andersen barostat的想法出发，给粒子的速度加了一个可以自由采样的scaling factor。这个因子作为额外的自由度引入到体系的Hamiltonian中，使得体系能量可以变化，但拓展后到pseudo-Hamiltonian 仍然是守恒的。

实际上各个热浴有各自的优点和缺点，下面我们通过一些例子来看不同热浴方法的特点。

加载依赖包。本文档使用的镜像为dmff_0.2.0-notebook。

Toy System -- Independent 3D Harmonic oscillator

Unit convention:

• Energy in $kJ/mol$
• Velocity in $ns/ps$
• Length in $\mathrm{ns}$
• Time in $\mathrm{ps}$
• Mass in $g/mol$ or $\mathrm{amu}$

体系描述 Description

让我们来假设体系有N个粒子。N个谐振子相互独立，体系的势能为： $$V(x,y,z)=\frac{1}{2}k(x^2+y^2+z^2)$$

$k$为弹簧弹性系数，假设为$100kJ/mol/nm$。假设体系与$200\mathrm{K}$的热浴耦合.

根据上述理论, $N=100$且$T=200\mathrm{K}$时体系的温度涨落应该为： $${\sigma }_T^2=\sqrt{T^2\ast \frac{2}{3N}}=16.3\mathrm{K}$$ It's a huge variance! 看看下面的实验是否与理论之相符合？

Stochastic Thermostat 1: Andersen Thermostat

在Andersen Thermostat中，我们需要随机抽取$m$个粒子让它们与热浴粒子相碰撞。因此算法需要执行以下步骤：

• 按照Velocity-Verlet算法演化位置和速度；
• 随机选取m个粒子（均匀随机，uniformly）；
• 改变m个粒子的速度，速度数值从Boltzmann分布中采样

这里面存在一个超参数，即粒子的耦合强度（碰撞频率 或 被选中的概率）。理论上，任意一个取值在$[0,1]$之间的碰撞概率都可以使温度稳定在期望值附近，并产生合理的系综分布。

但过于频繁的碰撞（即过强的耦合会导致所有的动力学信息损失）

Question: 一个粒子相邻两次碰撞的时间间隔符合什么分布？为什么?

（答：Possion Distribution）

Equilibrium is necessary!

现在我们得到了结果，让我们直接把结果画出来是什么样子的？

你会发现，系统的势能漂移了很久才达到我们想要的数值；而且似乎能量分布也不是很符合预期。

这是体系初始化造成的。

体系在初始化的时候，为了方便，起始点很容易出现在远离热平衡的地方，因此体系必须先 弛豫（Relaxation） 到平衡位置附近，我们才能收集到有效的数据点。

当然，理论上当采样时间无限长的时候，轨迹初始的一小段非平衡数据并不会对整体产生影响（因为此时他们的权重很低）。但由于我们只能做到短时间的有限采样，实际情况下这些非平衡数据存在过高的权重，会对数据分布产生很大影响。因此为了能够得到有效的数据，我们需要让体系先达到热平衡。

有两种方式可以选择：

• 预先模拟一段轨迹，不收集数据；
• 直接摒弃数据的非平衡部分。

其中，第二种方式在Monte Carlo模拟中也叫作 Burn in.

那么让我们看看扔掉前1/2数据后得到的结果吧：

现在基本符合预期了！！

注意，这里我们计算得到体系温度的涨落大概为$16\mathrm{K}$，这与理论值很接近。这也从数值实验中说明了，在NVT系综中，有限粒子体系在有限时间模拟下存在温度浮动。

Stochastic Thermostat 2: Langevin Thermostat

Langevin Thermostat在之前@hmj同学的notebook中已经提到，这里不再重复展开。

请有兴趣的同学移步：https://nb.bohrium.dp.tech/detail/2722437140

Darkside of Stochastic -- Missing of dynamics information

随机热浴的优势在于算法简单，能够很快实现系综的采样，没有复杂的数学推导，有直观的物理解释（粒子碰撞）。但它同样存在严重的缺点：无法准确描述动力学性质。

可以想象，粒子的轨迹应当是连续的。但是在Andersen热浴中，参与碰撞的粒子的历史速度信息被直接抹除，被assign了一个全新的速度。这种Markov Chain式的方式使得粒子全都在无规则的运动。在极端情况下（即每一步中所有粒子都参与碰撞），粒子不存在速度信息。

这对于一般的统计热力学物理量的采样没有影响，但是对于描述粒子运动的物理量有很大影响，比如：扩散系数 Diffusion Coefficient.

另外，由于Andersen热浴的速度/动量直接从Boltzmann分布中得到，体系整体的动量是不守恒的。

在某些时候，体系还会出现平动自由度过热，振动自由度过冷的情况，即Flying ice cube

Langevin Dynamics 保存了部分动力学性质，但是其收到随机种子的影响，会导致轨迹对随机种子产生依赖。

Angle from Diffusion Coefficient

According to Kinetics theory, the diffusion coefficient is proportional to mean square displacement (MSD):

$$⟨(x-\bar{x}{)}^2⟩=2Dt$$

现在我们来计算一下MSD.

因为下面的case用DMFF算的太慢，敬请期待OpenMM案例 :)

Global Thermostat 2: Nose-Hoover Thermostat

NH Dynamics 在之前@hmj同学的notebook中已经提到，这里不再重复展开。

请有兴趣的同学移步：https://nb.bohrium.dp.tech/detail/2722437140