各种系综下Monte Carlo粒子模拟的实现

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$为总粒子数。对于移动粒子,该步的接受概率与正则系综相同。若改变体积:

  1. 通常随机产生一个体积对数的增量$\ln V_n = \ln V_o+\Delta\ln V$,其中$\Delta\ln V$满足设定范围$[-c,+c]$上的均匀分布。
  2. 改变体积时不改变粒子的坐标分数,即体积改变是对系统长度尺度的整体缩放,所有的位矢也按照相同的长度比例缩放。

该步的接受概率为

假如势能$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在实际中用得很少。

参考资料

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