Monte Carlo方法中的Metropolis算法起始于对正则系综的模拟。Metropolis算法的核心是根据细致平衡(detailed balance)原则设置Monte Carlo循环中每步产生一个构型,该构型能否被接受的概率$acc(o\to n)$。这里$o$标记移动前构型,$n$标记移动后构型。
正则系综:NVT
每一步的接受概率是
在Markov Chain Monte Carlo方法中,通常只需要直接计算能量差值即可。例如一个粒子的移动使得该粒子与其他粒子的势能发生改变,但是其他粒子之间的势能没有改变。所以只要每一步计算该粒子在移动前后与所有粒子的势能的 差就可以得到 $U_n-U_o$。
恒温恒压系综:NPT
每一步由随机数决定改变体积还是移动粒子,改变体积的概率为$1/N$,$N$为总粒子数。对于移动粒子,该步的接受概率与正则系综相同。若改变体积:
- 通常随机产生一个体积对数的增量$\ln V_n = \ln V_o+\Delta\ln V$,其中$\Delta\ln V$满足设定范围$[-c,+c]$上的均匀分布。
- 改变体积时不改变粒子的坐标分数,即体积改变是对系统长度尺度的整体缩放,所有的位矢也按照相同的长度比例缩放。
该步的接受概率为
假如势能$U$与粒子对的距离有关,则缩放倍数与具体的距离的次方$-l$有关,容易推导 $U_n=\left(\frac{L_o}{L_n}\right)^l U_o$。注意若对对势进行了长程截断,则需要按照缩放后的截断距离重新计算能量补偿,因为能量补偿与距离的关系通常不是$-l$次方。总之,对于距离以不同次方依赖的势能项分别放缩即可,计算的开支比较小。
巨正则系综:μVT
每一步由随机数决定改变粒子数还是移动粒子,改变粒子数的概率为$1/N$,$N$为总粒子数的期望值,通常是一个预设的估值。对于移动粒子,该步的接受概率与正则系综相同。若改变粒子数,则随机决定增加还是减少一个粒子:
增加一个粒子的接受概率是
减少一个粒子的接受概率是
其中$\Lambda=\sqrt{\frac{h^2}{2\pi m k_B T}}$为热力学de Broglie波长是定值,$\mu$为化学势是给定值。
微正则系综:NVE
微正则系综Monte Carlo与以上各方法明显不同,区别在于不使用随机数决定一步的接受与否。当然每一步选取粒子和选择方向还是要使用随机数的,只不过每步移动要看系统的总动能够不够。
由于系统的总能量是固定值,从构型可以计算出总势能,总能量减去总势能就得到总动能。注意总动能不能为负值。所以每一步粒子移动,势能会发生改变,假如粒子移动使得总动能为负值,则不能移动。这里就完全没有用到随机数来决定是否接受这一步。
微正则系综Monte Carlo在实际中用得很少。
参考资料
- Frenkel, Daan, and Berend Smit. Understanding molecular simulation: from algorithms to applications. Vol. 1. Elsevier, 2001.