分子动力学模拟在原理上是描述微正则系综（NVE）下的动力学过程的方法。在正则系综（NVT）中，需要考虑一个外部的恒温环境，能够按照一定的方式改变系统中粒子的动能（速度）。可行的变更速度的思路主要有两种：

随机改变：将Monte Carlo过程加入到分子动力学演化中；
拓展Lagrangian：在系统的Lagrangian中加入参数，使得统计结果等价于NVT下的统计。

以下介绍这两种思路具有代表性的算法：Andersen Thermostat和Nose-Hoover Thermostat。这里主要关注其算法的实现。

Andersen Thermostat

算法

简单来说，Andersen Thermostat的实现就是在分子动力学的每一次循环（每一时间步）之后加入Monte Carlo过程：

对每一个粒子进行判定：有$\nu\Delta t$的概率改变速度。
这里$\nu$代表系统与环境（heat bath）的耦合强度，$\Delta t$为时间步长。显然$\nu$越大，改变的几率越大，Thermostat的效果越好；但是过大的数值会让分子动力学过程失去意义。
若判定某个粒子改变速度，则由指定温度下的Gauss速度分布随机给定一个速度。
这里指定温度下的高斯分布即 $N(\mu=0,\sigma^2=T)$，其中$T$为无量纲温度。

效果

Andersen Thermostat在适当的耦合强度$\nu$下或者足够长的时间之后能够给出合理的相轨迹，对状态物理量实现了恒温。当$\nu$太小而且时间不够长时，温度的改变不充分，会导致相轨迹离散化，不满足遍历性要求。但是这个不算大问题，只要合理设置$\nu$即可。

Andersen Thermostat不能正确给出过程物理量，如扩散速率$D$等，因为$\nu$的物理意义不够明确，所以存在任意性；而此任意性会明显地影响过程量。或者说，Andersen Thermostat因为将Monte Carlo引入分子动力学，所以失去了某些分子动力学的严格性。但是对状态量如压强等是严格的。

Nose-Hoover Thermostat

算法

Nose-Hoover Thermostat利用了能勢修一扩展的Lagrangian，引入一个额外的参数来实现恒温。粒子的运动方程修改如下

$\begin{align} \dot{\brm{r}}_i&=\frac{\brm{p}_i}{m_i}\newline \dot{\brm{p}}_i&=-\frac{\pd U(\brm{r}^N)}{\pd \brm{r}_i}-\xi\brm{p}_i\newline \dot{\xi}&=\frac{1}{Q}\left(\sum_i \frac{\brm{p}_i^2}{m_i}-\frac{L}{\beta}\right)\newline \dot{s}/s&=\xi\label{eqn:1} \end{align}$

这里$U(\br^N)$是某时刻系统的总势能，$\xi$和$s$为引入的参数，$L$在三维$N$粒子系统中取$3N$；$Q$类比一个质量参量，代表系统与环境的耦合强度。其中式$\eqref{eqn:1}$是多余的，因为前三式已经完备。通常式$\eqref{eqn:1}$用来检验代码：三维$N$粒子系统的Nose-Hoover动力学保证$H_\text{Nose}$是守恒量

$H_\text{Nose}=\sum_{i=1}^N\frac{\brm{p}_i^2}{2m_i}+U(\br^N)+\frac{\xi^2Q}2+3N\frac{\ln s}{\beta}$

效果

与Andersen Thermostat相比，Nose-Hoover Thermostat最大的优势在于能够得到正确的过程量，如恒温扩散系数等。耦合强度$Q$对这些量影响很小。

Nose-Hoover Thermostat的致命缺陷在于，当系统受外力（即外势能）时，Nose-Hoover Thermostat不能给出正确的相轨迹，严重地违背了遍历性，最终给出错误的结果（在这种情况下Andersen Thermostat仍然适用）。这个问题的解决方案是使用一系列耦合的Nose-Hoover Thermostat来满足遍历性。

Nose-Hoover Chain

算法

在Nose-Hoover Thermostat的基础上，增加参数，相当于增加多个耦合的Nose-Hoover Thermostat：

$\begin{align} \dot{\brm{r}}_i&=\frac{\bm{p}_i}{m_i}\newline \dot{\brm{p}}_i&=-\frac{\pd U(\brm{r}^N)}{\pd \brm{r}_i}-\frac{\eta_1}{Q_1}\brm{p}_i\newline \dot{\xi}_k&=\frac{\eta_k}{Q_k},\qquad k=1,...,M\newline \dot{\eta}_1&=\left(\sum_i \frac{\brm{p}_i^2}{m_i}-\frac{L}{\beta}\right)-\frac{\eta_2}{Q_2}\eta_1\newline \dot{\eta}_k&=\left(\frac{\eta _ {k-1}^2}{Q _ {k-1}}-\frac{1}{\beta}\right)-\frac{\eta _ {k+1}}{Q _ {k+1}}\eta_k,\qquad k=2,...,M-1\newline \dot{\eta}_M&=\left(\frac{\eta _ {M-1}^2}{Q _ {M-1}}-\frac{1}{\beta}\right) \end{align}$

这里$\eta_k$类比动量；$Q_k$类比对应的质量参量，代表各个Thermostat之间的耦合强度。三维$N$粒子系统的Nose-Hoover Chain动力学保证$H_\text{NHC}$是守恒量

$\begin{equation} H_\text{NHC}=H(\br,\brm{p})+\sum_{k=1}^M\frac{\eta^2_k}{2Q_k}+3N\frac{\xi_1}{\beta}+\sum_{k=2}^M\frac{\xi_k}{\beta} \end{equation}$

效果

通常在使用两到三个耦合的Nose-Hoover Thermostat，就可以得到很好的结果，可以正确处理单个Nose-Hoover Thermostat不能处理的外势能体系的问题。相应地，计算量也会随着耦合数量增加而增加。

参考资料

Frenkel, Daan, and Berend Smit. Understanding molecular simulation: from algorithms to applications. Vol. 1. Elsevier, 2001.