Metropolis 蒙特卡罗方法、动力学蒙特卡罗方法、分子动力学方法这三种模拟方法有何特点与差异？ 第1页

偶然发现这个问题以及诸位大神的回答，大致看了一遍之后发现诸位的回答虽然没啥大问题，然而无论是MC还是MD，本质上都是一门技术而不是一种理论，只要是技术就必须勤加练习才能真正掌握，光说不练看再多书也没用。而上面的大多数回答都太过理论，显然对于刚入门的同学是很难看懂的。而且大部分回答也过于罗嗦，毕竟问题问的是三者的“特点”和“差异”，如果只是分别论述三者的特点，其实等于什么都没回答。

从技术的角度来讲，想要真正理解这个问题，要只要明白MD（Molecular Dynamics分子动力学）和MC（Monte Carlo蒙特卡罗）的最大差别即可：MC只需要算势能（Potential Energy），MD除了算势能还需要算力。就是因为这一最大差异，导致MC和MD的应用范围不同，特点也不同。

因为可以算力就可以算速度，所以MD的可以由当前的坐标和速度计算出下一步的坐标和速度，也就是说MD可以得到时间连续的轨迹。

因为没有力也就没有速度，所以MC只能靠随机移动来决定下一步的位置，然后通过势能的不同计算每一步的权重，或者根据前后两步势能差的大小来决定是否接受这一随机移动的概率（Metropolis MC)。

从统计学上来讲，两种方法对体系的抽样是等价的。然而MD相比MC可以得到时间连续的轨迹，所以MD可以计算的性质更多。而且MD是根据力按部就班地一步一步计算，即使是对于复杂的真实体系也很容易实现MD软件的通用化。相比而言，因为复杂体系的随机移动是非常困难的（越复杂的体系在坐标随机移动之后就越容易出现重叠，这一点是随机移动所不允许的），基本上不同的体系需要完全不同的随机移动算法，所以MC软件非常难以实现通用化。

然而MD必须计算力这一点同样也制约了MD在某些领域的应用，因为力的计算需要势能对体系坐标是连续可导的，然而并不是所有体系的都可以满足这个条件。而MC则不需要计算力，而且因为MC的实现远比MD简单，虽然没有通用的软件，针对特定体系编写一个MC程序往往也不复杂。所以MC特别适合自建模型体系的模拟，因为模型体系的势能有可能十分复杂，无法或者很难求导，而且模型体系的随机移动算法也更容易实现。

综上所述，MC和MD最根本差别就是是否算力。需要算力的MD程序更加复杂，但是可以实现通用化，且可以得到时间连续的轨迹，所以适合于模拟复杂的真实体系。而不需要算力的MC程序则比较简单，即使体系势能无法求导仍然可以模拟，但是随机移动的算法难以实现通用性，导致MC程序也很难实现通用性，所以MC更适用于模型体系。

至于KMC（动力学蒙特卡罗）我本人并没有真正用过，所以无法回答的太深。然而KMC本质上仍然是MC，使用MC计算不同“态”之间的势垒从而间接获得“动力学”，然而却仍然无法得到时间连续的轨迹。

这个问题很大，详细展开足够在一本分子模拟的教材里写满几章了，了解它们最好的方法应该是去啃书，比如分子模拟的经典教科书："Understanding Molecular Simulation, From Algorithms to Application" by Frenkel & Smit。

这里先讨论MD和Metropolis MC, 暂不讨论Kinetic Monte Carlo (KMC)。这是因为传统的MD和MC更具可比性（为简单起见，这里MC就是指狭义上的使用Metropolis算法的蒙特卡罗，而这里谈到的MD仅限于平衡态的MD模拟，暂不讨论非平衡态的NEMD），两种方法产生的初衷都是为了计算统计力学里的系综平均：假设系综所对应的相空间的概率分布为（p为动量，q为坐标），那么对于任意平衡态热力学量M, 它的系综平均 可以表示为

(1)

说得更直白一点，传统的MD和MC基本上是在干同一件事：算积分。只不过最早的MD是在微正则系综（NVE系综）里算积分，而最早的MC是在正则系综（NVT系综）里算积分。

先说MC。众所周知，正则系综里相空间的概率分布就是玻尔兹曼分布

(2)

其中是体系的哈密顿量， 是配分函数， . 以一维线性谐振子为例，这个相空间概率分布就是下图所示，其实就是个高斯分布。

知道了概率分布，那么正则系综中任意平衡态热力学量M的系综平均 就是下面这个积分

(3)

对于（3）式中的积分，绝大多数情况下是没有解析解的，因此只能通过数值求解。最简单粗暴的用蒙特卡洛方法求解（3）式积分的方法就是在整个积分区域上均匀地随机打点，求得每一点处的值，然后对所有点求和再除以配分函数即可。例如对于一维线性谐振子，我们可以如下在相空间上均匀地打点

这种算法，的确非常简单，但效率很低。因为位于相空间中的值非常小的区域，系统处于这个区域的概率非常低，也就是说这些区域大体上来说对系综平均 的贡献是非常小的，然而在采样过程中，由于我们是随机均匀打点，这样的值非常小的点和的值非常大的点的权重是一样的。换句话说，我们浪费了不少时间在那些权重很小的点上。那么是否有这样算法，能够使产生的点的分布符合玻尔兹曼分布，也就是说，在的值非常小的区域少打点，在的值非常大的区域多打点（如下图所示），然后把每一点处的M值相加，就能得到我们需要的系综平均呢？Metropolis算法就是为了实现这样的目标产生的。关于Metropolis算法的细节，我不在这里展开，很多书籍里都有描述。

再说MD。MD是基于牛顿方程的，而从牛顿方程是可以很自然导出孤立体系在保守力作用下是能量守恒的，因此最初的MD是在微正则系综(NVE)里计算系综平均。根据统计力学的等概率假设，对于总能量为E的微正则系综，其相空间概率密度分布为

(4)

其中是体系处于能量为E时的态密度(density of state)

(5)

那么热力学量M在微正则系综里的系综平均 就是

（6）

其中V是积分的区域，也就是当能量为E时系统所占据的相空间。拿一维线性谐振子作为例子，当能量为E时其在相空间中区域对应于下图中的圆环，那么这个系综平均就是根据（6）式在对应于能量为E的圆环区域上进行积分。MD所做的，就是让粒子跑遍这个圆环上的每一个点（历经各态，ergodicity)，把每一点上M的值加起来求平均，得到的就是系综平均(time average equals ensemble average, this is what we call ergodicity)。

综上，MD和MC都可看作是为计算系综平均的采样方法。采用Metropolis算法的MC不局限于在正则系综里的采样，可以应用于其他系综，比如NPT和巨正则系综。至于微正则系综，理论上也可以使用Metropolis算法，只是由于各个态之间能量相等，概率也相等，这时使用Metropolis算法不能起到加速采样的作用，等于白用。同样地，MD的采样也不局限于微正则系综，当把一个孤立体系和一个恒温器(thermostat)或者恒压器(barostat)耦合起来，就可以进行正则系综或者NPT系综里的模拟。当然，如果耦合的恒温器或者恒压器是以extended Lagrangian形式实现的话（比如Nose-Hoover Thermostat), 广义上来讲MD模拟还是可看作在一个大的微正则系综中进行。

MD和MC的比较

1. 算法的难易程度。一般说来，MD的算法比MC要复杂一些。这是因为，对于MC而言，每产生一帧新的构型只要计算能量然后根据acceptance ratio移动粒子就足够了。但对于MD，每产生一帧新的构型所需要的不只是能量，还需要每个粒子的速度和受力来对牛顿运动方程进行数值积分，这就加大了计算量。如果还要考虑为了提高计算效率的Verlet neighbor list等设置，MD的算法会更复杂一些。
2. 所能计算性质的多少。很明显，MD能比MC计算更多的性质，MC只能计算一些与动力学无关的平衡态性质，比如能量、比热容等等，MD既可以计算与动力学无关的性质，也可以计算与动力学相关的性质，比如扩散系数、热导率、黏度、态与态之间的跃迁速率，光谱等等。这是因为，使用Metropolis算法的MC，为了快速得到平衡态的分布，把整个系统简化为一条马尔可夫链，第N个构型只和第N－1个构型有关联，与第N－2或者更早的构型没有什么关联。在实际情况中不是这样的，系统是有记忆的，在时刻处于相空间中点的系统与在时刻处于的系统是有关联的，它们之间的关联可以用条件概率来描述。很多动力学性质都和这样的条件概率有关系，而MD产生的数据是可以计算这些概率分布的。所以世界上没有免费的午餐，MC虽然算法简单，但大体上只能计算静态性质，如果要计算动力学性质还是要上MD。
3. 可使用的软件。可使用的MD的软件远多于MC的软件。我个人觉得很大程度上是由于MD能够计算的性质远比MC多，并且数据可视化后给人更多直观的感受。比如说，你可以把MD模拟蛋白质折叠的轨迹文件用软件可视化之后做成动画，动画就比较直观地演示了蛋白质结构随时间的变化。可是如果你用MC做同样一个体系的模拟，你把轨迹可视化之后看到的可能就是一堆原子像无头苍蝇一样乱动，最后稳定在某个构型附近。一般来讲，MD软件的开发也不难，一旦把积分牛顿方程的engine写好，往里面添加新的相互作用势并不难，这样每添加一点新的功能，MD软件就可以适用于更多的体系。相比起MD软件的百花齐放，分子模拟中学术界通用的MC软件并不多，很多课题组都是自己写MC代码，但适用范围仅限于自己研究的体系。
4. 增强抽样(enhanced sampling)。这个问题 @JimKarrey 在他的回答中已经提到了。 这是分子模拟中的一个大坑，不管MD还是MC都会面临在系综中采样不足的问题。这方面知乎上讨论已经很多，可以参考下面的链接，在此不做详细描述。

最后说一下KMC，我在这方面不是很懂，期盼有高人来解答。

Kinetic Monte Carlo(KMC)和主要模拟平衡态性质的MD和Metropolis MC相差就比较远了，它主要模拟系统随时间的演化而非平衡态性质。KMC的基本思想是把整个系统划分成若干个态(state)，并获得态与态之间单位时间的跃迁几率（这个可以通过分子模拟或者过渡态理论来计算），然后通过具体的算法决定每一步往哪一个态跃迁，就可以模拟系统在各个态之间随时间演化。我没有做过KMC的研究，并不太了解KMC的算法，仅仅粗略看过一两篇用KMC模拟扩散的文章。以下的例子是我粗浅的理解，可能不太准确。

我们想用KMC模拟粒子在二维晶格上的扩散，并求出扩散系数。假定在单位时间步长内，粒子只能跳到它最相邻的格点。如图所示，粒子某时刻处于A格点，那么粒子下一时刻只有可能跳到B,C,D,E格点，单位时间内粒子从A跳到这些格点的跃迁几率分别为。结合这些跃迁几率，KMC的算法就可以通过生成随机数来决定粒子下一步所在的格点。在模拟足够长时间后我们可以算出粒子的均方位移(mean square displacement, MSD)， 均方位移和扩散系数成正比关系，因而可以求出扩散系数。

当粒子只能在一维情况下运动（即只能从A到B或者从A到D), 并且， 这种情况下KMC其实就是在模拟一维随机游走(random walk)。

KMC的具体应用应该还很多，比如晶体生长什么的。

这个问题非常好！我之前也想过一段时间，下面就把我的理解说一下。

还是先说一下为什么我觉得这个问题好。

1. 在利用统计物理的方法求系综平均时，肯定需要采样，而采样方法很多，到底应该用哪一种，这些方法的联系在哪，是否都是一种方法在不同条件下的应用（实际上都是MCMC）。
2. 上学期学了一些统计学习，里面也用到了很多采样方法，实际上，如果能把统计物理里的采样方法的理解提升到更高层次，完全是可以理解很多机器学习的问题的。传统机器学习里的Kalman滤波实际上就是用了Langevin动力学（这里的分子动力学也是一样的），Particle滤波用到的就是在重要采样的基础上发展起来的DMC。特别火的AlphaGo基于的RL，实际上就是对不同形式的Bellman方程进行采样，对价值函数进行估计，然后采用不同的策略，比如贪心算法， -贪心算法。怎么估计价值函数、做出决策，又发展了一系列的，比如Deep Q Network，Monte Carlo Tree Search。统计物理里在样本较少情况下常用的bootstrap方法也在RL中有体现，比如experience reply，这些都是AlphaGo算法里很核心的部分。

要理解这三个算法的不同，需要对概统里的采样方法有一点了解。一般的采样方法有

1. 直接采样：把累计分布函数的反函数拿来采样，但是C.D.F.一般不太容易得到，所以不常用。 ，通过均匀分布采样
2. 接受-拒绝采样：拿一个大的函数把待采样的分布套起来（因为有归一化要求，必须要用一个大于1的数乘上某个分布） ， s.t. 。采样按照拒绝率判断是否接受采样样本。采样位置 ，接下来判断是否接受， ，如果 ，接受x，否则，拒绝x。相对于直接采样法效果好一些，但是，要找到好的G(x)不是简单的事情。如果c太小了，很难满足要求，如果c太大了，接受率太低，采样效率不高
3. 重要采样：通过变换 ，通过按照g(x)采样，对带权的函数求平均来求之前p(x)下的期望。重要采样是为了解决p(x)采样的样本点比较密集处f(x)的值较小导致的偏差较大的情况。但也有不足，就是g(x)不好找，并且，在某些情况下，采样次数要比较多才能收敛。
4. MCMC的方法：这里的三种方法都属于MCMC，Markov Chain Monte Carlo。都是要构造一条Markov链，使其稳定分布收敛到待采样的分布，这样，收敛之后每一步的结果都是从目标分布p(x)中采样的样本。Markov链能收敛的必要条件就是细致平衡条件成立，即

下面就分别说下题主提到的三种采样方法。

用到的一些关于微观动力学方程的知识可以在专栏里找到。

一、Metropolis

1. Metropolis-Hasting算法：Metropolis算法实际上是Metropolis-Hasting算法的特殊情况。MH算法的核心在于，我们通过提议分布来控制采样位置，接受率来控制是否进行采样，跟接受-拒绝采样比较像。假设有便于采样的提议分布 。在当前采样时刻t，需要知道上一步t-1的样本点x，t时刻的样本点 。均匀分布随机数 与接受率 进行比较，如果 ，接受y，否则，拒绝y。
2. Metropolis算法：当提议分布是对称的情况，i.e. ，有
3. 物理中的Metropolis算法：由于在物理里，我们的提议分布都采用高斯分布等对称分布的形式。所以，MH算法退化到Metropolis算法。
4. 应用：用Metropolis算法计算Ising模型，每次是否接受自旋反转，都是通过稳定分布（正则分布）之比来验证。

二、KMC

1. 为什么要引入KMC，Kinetic Monte Carlo。回忆一下之前说的MH算法相当于接受-拒绝采样，而接受-拒绝采样的不足之处在于接受率会不高，比如当模拟退火时，温度太低，可能会导致接受率非常低。
2. KMC的优势：实际上KMC相当于Metropolis算法在时间上进行粗粒化。每次不判断是否接受提议，而是通过轮盘赌的方法选出一个一定接受的操作，并且计算这个操作所需要的时间。相当于，在Ising模型的例子中，有可能1000步才会接受一次自旋反转，但是在KMC里不跑1000步，1步接受，并且估计这一次操作所需要的时间为1000步。
3. 应用：用的比较多的是在化学反应系统中。系统中有 个反应通道，每个R对应一个化学反应方程，r是该反应的速率，根据浓度和速率常数计算。轮盘赌抽样 返回满足条件的k，则这一步，反应 发生，其对应的分子相应变化，这一步的时间，

三、分子动力学, Molecular Dynamics, MD

1. 首先要解释一下，题主在提问里说，MD是确定性算法，这是错误的。首先想一下MD怎么工作的，是不是用到了高斯随机数。其次，MD会被误认为是确定性算法，还是由于Langevin方程非常类似于牛顿方程，这一方面是一件好事，即是你不懂Ito公式，不懂随机过程，也可以通过牛顿力学来理解MD，创造了不少科研岗位，但另一方面，门槛低了，有很多人误解做计算，特别是MD不需要什么数学物理知识，只是用用软件，调调参数。
2. Gibbs算法的引入：首先介绍另一个MCMC算法，Gibbs算法，不是Gibbs发明的，是在研究Gibbs随机场的情况下发明的。之前说过MH算法的不足在于有接受率，所以需要改进，引入了KMC，在时间尺度进行粗粒化，但是从上面的例子也可以看出，KMC能采样的系统是不够精细的，上面例子中你并不知道到底哪个分子在什么位置发生了化学反应，只知道系统的浓度发生了变化。因此，在应用中，KMC是不够的。另外，由于分子系统的自由度非常大，相当于要在6N维的相空间上进行采样，接受率即使是0.99，也会产生维数灾难，使得整个样本向量的采样接受率几乎为0。要解决接受率的问题，就要从接受率本身来看，如果能使提议分布满足 ，接受率就是1。Gibbs算法就是从这个角度出发。特别的，如果你的提议分布就是Markov链的转移概率，因为MCMC中的Markov链满足细致平衡条件，Gibbs算法就成立了。
3. Gibbs算法：记待采样的样本向量为 ，记t时刻采样时，除去第j个分量的向量为 ，采样过程中，在t时刻， ， 为t时刻采样得到的样本向量。这个算法实际上就是在 中沿各个坐标轴分别随机游走。
4. 应用：在MD中，提议分布采用高斯分布，因为从Langevin方程离散得到的向前差分采样就是从高斯分布里采样。稳定分布则通过Langevin方程对应的Fokker-Planck方程得到。这两个概率密度是满足细致平衡条件，也就是Gibbs算法中对提议分布的要求，所以，每次采样都接受。因此，在MD中每一步构型改变都是在相空间里进行了一次Gibbs采样。

四、补充点东西

1. 前面提到的几种基本采样方法，实际上在统计物理里也有很多应用，比如重要采样就衍生出了像DMC (Diffusion MC)等方法，在高维空间采样中有比较重要的应用。不过用DMC还是有很多要注意的技巧的，如果不注意，DMC的采样结果非常容易收敛到退化分布。
2. 三种方法的特点：Metropolis在Ising模型及其类似的体系里用的很多，好处是微观细节反映的比较好，但是因为接受率的问题，采样速度慢，并且只在低维空间比较有效。KMC在有平均场近似的地方用的很多，比如复杂系统、反应扩散系统、用统计物理研究社会问题。采样速度比较快，在要研究的时间尺度比较大时可以考虑，但是因为只体现了平均场近似之后的性质，对微观过程的体现不太好，而且也只适用于低维空间。MD对微观细节的体现非常好，适合研究像polymer、protein这样的生物物理模型，因为没有拒绝率的问题，在高维空间采样的效率和稳定性都比较好。并且，MD能最大化的得到系统的几乎所有力学量，因为这个采样方法几乎是在系综（统计物理里对概率空间 的描述）上给出了样本的足够多次实现（Markov链收敛之后，想要采样多少次都可以），根据强大数律，力学量的样本均值a.s.收敛到期望。
3. MD的基础上有很多发展，比如各种保辛算法的应用，通过Markov链轨道概率的分解发展出来的Path Integral MD，在时间、空间上的粗粒化方法。

被李老师 @Yongle Li 指出微正则系综里的分子动力学是确定性算法，i.e. 哈密顿力学。不过个人理解，实际上是Langevin方程的特殊形式。由于微正则系综下不考虑能量的涨落，所以扩散过程因扩散张量 而退化到确定性过程，而扩散张量也是噪声的关联张量，因为其为0，噪声也就退化为取值为0的退化分布。因此就不需要考虑噪声项。Langevin方程退化为 ，在确定性过程中，根据Liouville定理，概率密度对应的Fokker-Planck方程的解为退化分布，取值在前面常微分方程的解上。而前面的方程即为哈密顿正则方程

部分内容见专栏