中国矿业大学分子模拟课件第十二章.pdf
分 子 模 拟 牛继南 njn0516cumt.edu.cn 2011.4 第十二章 蒙特卡罗计算方法 12.1 基本原理 广义定义：指凡是在计算中引用随机数的 计算法，英文名称Monte Carlo，简称MC。 以计算半径为单位长度的圆的1/4面积为 例说明蒙特卡罗方法基本原理。 选取一组随机点P(x,y)，它们的范围为 0x1，0y1。 如果此点满足 则被接受，否则被排除。 如果随机点数为N，被接受的点数为n，则 我们利用随机点求得的1/4圆面积为： P点落入1/4圆面积的概率乘以1×1正方形面 积。 1/4圆的真实面积为·r2/4=0.785398，而不 同随机点数得到的面积为： 以上的问题相当于用蒙特卡罗法求积分： 可以将这种方法推广求解多重积分，一般 用于普通方法难以求解的积分，如对于： 计算时，在变量xi可能的区间内任意取n组 随机点 ,i=1,n. 相当于计 算： 这种方法的误差为n-1/2，故选取的点越多， 计算的积分值越精确。 对于积分函数： 可以改写为： 其中p(x)为配分函数，设在范围x1,x2内按 照p(x)分布选取n个随机数x，则依照统计学 积分式可写为： 式中表示对n个随机数的平均。 对于正则系综中任意坐标，函数F的平均值 为： 式中P(x)为： 为系统位于特殊结构 的概率。 定义结构积分为： 则F的平均值为： 理论上若任意选取大量的结构点，则可按 MC法求解积分，一般步骤为： 1. 随机产生N个原子坐标； 2. 计算此构型的势能和函数F； 3. 计算玻尔兹曼因子 4. 累加玻尔兹曼因子，并将函数与玻尔兹曼因子相 乘后累加，回复步骤1； 5. 当计算Nt循环后，求出F的平均值为： 理论上可以用蒙特卡罗方法计算统计的平 均值，但若计算大体系必须选取大量的样 品点，而这些样品点所对应的结构多数因 为能量高，因此不容易出现。此外，依靠 此方法计算积分时，必须先计算结构积分 ，但实际上此积分不易收敛。 基于这些原因，事实上直接用MC法并不实 际，通常采用Metropolis的重要取样法。 12.2 重要取样法 以NVT系综的某一物理量A的平均值为例， 统计力学算式为： 式中NVT为概率， 为对应的结构。从统计 角度讲，有些结构出现的概率大（trial）， 有些结构出现的概率小，而重要取样法计 算平均时设法选取重要的 ，使积分仅为 这些重要贡献的求和： n为重要取样点。 这种方法是Metropolis，Rosenbluth，Teller 等在1953年首先提出，被称为Metropolis法 ，是现在MC计算的基础。 此法的重点是建立系统的马科夫链(markov chain)，使最终满足重要的取样点的要求. 马科夫链符合两个条件： 1. 设系统个状态所成的有限集合为 则每次选择的随机数取样点产生的状态 必为此 集合的一个状态。 2. 每次所选的取样点仅与前一个取样点有关，而与 由何处开始取样无关。 以计算机运行为例： 若当天正常运转，则次日正常运转的概率 为60%；若当天停机，则次日也停机的概率 为70%。 考虑此问题时分为正常（状态1）和停机（ 状态2）两种状态，由一种状态转换成另一 种状态的概率用ij来表示，则： 11表示第一天正常第二天也正常，为0.6 12表示第一天正常第二天停机，为0.4 21表示第一天停机第二天正常，为0.3 22表示第一天停机第二天也停机，为0.7 明显有： 可以建立转换矩阵 设1为系统处于状态1的概率，2为系统处 于状态2的概率，则系统的状态概率可用概 率向量表示： 设 为系统起始状态的概率向量， 为系 统经过第一次转换后的概率向量，即: 再经过一次转换后: 以此类推，即得到马科夫链。 如果第一天计算机正常运转和停机的概率 相等，即 =（0.5,0.5），则第二天的概率 向量为： 即第二天正常概 率为0.45，停机 概率为0.55 第三天的概率向量为： 经过多天的计算得到极限值： 若概率达到极限，则再乘转换矩阵必得到 相同的概率向量，如： 这时表明系统的概率分配值不再变化，系 统达到平衡。 推广到N中可能的系统，则 为N维向量， 转换矩阵 为N×N矩阵，矩阵中每行元素 的和为1，即 满足极限概率值的条件为： 为了得到这个方程的解，通常采用微观可 逆条件来替换，即： 如果对m的所有态求和，则 = n m nmnnm m nmn m m 即为上面的式子。 第一个提出解这个方程的为Metropolis等， 他们引入对称矩阵，即mn=nm，表示选择 mn和选择nm的概率相同，设 分 别为状态m和状态n被接受的极限概率， 若已达到极限状态，且 ，则m转换为n 的概率等于选择mn的概率，即 若 则m转换为n的概率等于选择 mn的概率乘以n状态与m状态被接受的概 率比。 同时还需满足同一状态的转换概率： 由于多粒子系统有不同的状态，故转换矩 阵十分庞大，但极限状态可以由统计力学 得到。对任一结构 ，其极限概率为 因此 式中， 为二状态n与m的势能差。 当n状态的势能低于m状态的势能，则此结 构的转换可以被接受，若n状态的势能高于 m状态的势能，则根据前面的式子来接受此 移动。 实际当中可将玻尔兹曼因子与一个在0-1之 间的随机数比较，如果大于次数则被接受 ，如果小于次数则排除此转换。 按照这个图，随机数在左边靠近坐标为0的 时候，配分函数值较大，容易出现在函数 下方，所以容易被接受（此时虽然UnUm ，但能差很小）；如果随机数在右边，配 分函数值接近于0，因此随机点易出现在函 数上方，被排除（此时UnUm，且能差很 大）。 因此通过这样的转换可排除不可能的移动 ，达到重要取样的作用。 转换过程中转换矩阵 随系统的不同而不 同。 12.3 Metropolis具体计算步骤 以二维空间为例，系统中有数 个粒子，处于系统m。产生另一 状态时随机选择系统中的一个粒 子i，将其由原来的位置 以相同 的概率移动到一定范围R内的任一位置 。 范围R的中心为 ，边长为 。这样可以 产生NR个新位置，即NR个新的状态，因此 选择mn的概率为1/NR，可定义 为： 下一个步骤为计算移动的玻尔兹 曼因子。设选择的粒子i由原来的 移动至新的位置 ，设系统中 的粒子作用仅为LJ势，则： 即只需计算能量变化量。若 ，将玻 尔兹曼因子和随机数 比较，若大于 次数则被接受，反之则被拒绝。若接受以 新的位置为i原子位置，否则回到原位。 对含有N个粒子的NVT系统，Metropolis的 计算步骤如下： 1.任一挑选一个粒子，产生任意方向的位移 ，位移的大小为一定范围内的任意值，设 该粒子坐标为 ，移动后的坐标 为 则： 式子中Rand-1,1为随机数， 为计算中的最 大位移。 2. 设移动前系统的势能为 ，移动后系统的 势能为 ，则 若U0则此移动完全被接受，被接受的概 率为1，若U0则被接受的概率为 当U0是，可另外产生一个范围0-1之间的 随机数p，若 则此移动可接受 概率为 ，如果 则此移动无效，粒子坐标不变。 重复1和2的步骤，收集被接受的结构，并 与每一接受的结构计算一些热力学或结构 的物理量的平均值，如总能量、均方根位 移、次序参数等。如果这些量收敛，则系 统达到平衡。利用这种方法可将许多不重 要的结构省略，增加计算效率。 通常在计算时调整最大移动范围的值，使 所产生的移动约有50%可被接受。 传统的Metropolis计算方法每次仅选择一个 粒子移动。 Chapman提出32粒子同时移动的MC计算。 12.4 各类系综的蒙地卡罗计算 1. NPT系综 MC方法可以依据实际的需要，用于计算各 种不同的系综。NPT系综的体积可以在计算 中改变，达到平衡时系统的平均密度趋于 定值。设NPT系综为含有N个粒子的立方体 盒子，则任何物理量的平均值为： 式中ZNPT为单位立方体的结构积分，V为系统的体积， 为一组度量化坐标，其 中 , 为系统的边长。 系统的任一结构 其极限概率按统计力学 可求得： 以Metropolis方法产生(NPT)系统的马科夫链 的计算步骤为： (1) 随意移动一个粒子的坐标，或随意改变系 统的体积，或两种步骤都进行，即 为1-0之间的随机数， 为此类随机数组成的项链 ， 。 控 制向量化坐标和盒子的最大该变量。一般计算中选择的 使新位置的接受率为30%-50%。 (2) 产生新的坐标后，计算系统焓的改变量 即： 若 ，则移动完全被接受，若 则比较 和0-1间的随机数p，以 决定此移动是否被接受。 重复（1）和（2）步骤，并计算系统体积、焓的 平均值，以检验其是否收敛。由体积、焓的统计 涨落可计算一些相关性质。例如利用系统的位力 可以计算压力，此计算的压力和设定的压力是否 相等，可以用来检验MC计算的正确性。 2. 巨正则系综 巨正则系综 为体积、温度和化学势 皆为定值的系统集合，系统的粒子数为可 变化的参数，任意物理量A的平均值为： 式中， ，为系统的活性， 为德布罗易波长，h为普朗克常数， 为一 组度量化的坐标。 配分函数为： 按照此式系统物理量的平均值与粒子数目 有关，但此粒子数并不是定值，因此求物 理量A的平均值公式中为加和号而不是积分 号。 对于系统任意结构的极限概率为： 这类系统通常采用Norman和Filinov所发展 的方法。 产生马科夫链的计算步骤为： （1）任意移动系统中的一个粒子， 并依照一般的Metropolis方法决定此移动是 否被接受。 （2）任意消除系统的一个粒子，消除前后系 统的极限概率比为： 转换为活性 消除函数 若消除函数 则此消除步骤可完全接 受，如果 则比较 和随机 数的大小来决定此粒子的移动或体积的改 变是否被接受。 （3）随意产生一个新的粒子，产生前后系统 的极限概率比为： 这里定义 为新生函数，由此函数和 Metropolis方法判断此步骤是否被接受。 重复（1）（2）（3）的步骤，并计算系统 的构型能量、压力及密度以判断计算的收 敛性。 这种方法可以直接计算系统的自由能： 上面的消除和产生步骤，也可以由以下的 比较式来决定： 12.5 进阶蒙地卡罗计算法 一般的MC计算中，每次仅随机选取一个粒 子随机移动，因此系统中所有粒子的移动 概率都相等，但这种移动方式未必是最有 效的移动方式，例如经常有分子中的某些 粒子移动概率大于其它部分，或是移动粒 子时最好能使得粒子的移动遵循某些方向 。 1. 优先取样蒙特卡罗计算法 由Owicki与Scherage首先提出，为一种改进 取样点的计算方法。 以稀薄离子水溶液为例，此系统中最重要 的作用为溶质与溶剂间的作用。溶剂间的 作用不重要，溶质间彼此的作用更低。因 此考虑溶质与溶剂的作用，则以第一媒合 层的溶剂分子与溶质分子的作用最为重要 ，在此以外的溶剂分子，相距较远重要性 降低。 按此图，定义围绕溶质分子的特定区间为 Rsol，以此区分此区间内外的溶剂分子，则 可以定义一个介于0与1间的参数p来决定分 子的移动频率。 计算步骤如下： （1）在系统中随机选择以溶剂分子。 （2）若选取的分子处于Rsol区间内，则可以 产生移动。 （3）若所选的分子处于Rsol区间外，则产生 一个介于0到1之间的随机数，如果所选定 的p值大于或等于所产生的随机数，则移动 所选定的分子。否则，分子不移动。 （4）若分子不移动，则回复到步骤（1） 2. 力场导向蒙地卡罗计算法 此计算法是Pangali等提出的，此方法是以 粒子在系统中所受的力的方向为移动粒子 的依据。一般的Metropolis计算都是随机决 定所选粒子的移动方向，因此往往造成所 选的结构能量过高而不能被接受，以至效 率很低。对于含N个粒子的LJ系统为例，设 粒子i和j间的作用势为ij，则力场导向计算 法的步骤为： （1）随机选择移动的原子 （2）计算此粒子所受的力及方向 （3）由步骤（2）中力的方向决定所选粒子 的移动方向 式中Rand-1,1为随机数，rmax为计算设定 的最大移动距离。 力场导向计算法对氢键作用的分子系统特 别有效，因为此类系统中分子间具有一定 方向性的作用力，用引力场导向计算法可 以节省许多取样的浪费。Rao等将此方法应 用于计算纯水系统，并证实力场导向MC计 算结果较一般的MC方法精确，与分子动力 学计算的结果接近。 3. 明智的蒙特卡罗计算法 MC计算是基于随机取样与随机移动的原理 。前面的力场导向MC计算中合理的决定了 粒子移动的方向，移动的幅度有规定的最 大移动距离和随机数决定。如果将系统得 到一些区域特性一并考虑，则可以计算结 果的精确性。如计算溶液时，将溶质分子 的运动考虑为分子自身力场所产生的运动 和由溶剂影响而产生的布朗运动。这种方 法被称为明智的蒙特卡罗方法，由Rossky等 提出。 明智MC和力场导向MC类似，计算中都考 虑到粒子的受力及整体环境的影响。这两 种方法的计算结果远优于传统的MC，但计 算的复杂程度却几乎与MD和BD接近，所 需的时间也接近。 由于MC计算只能获得统计的平均值，而无 法获得系统动态的资料，因此通常情况下 用的并不是很多。 MC参考书： 1.Monte Carlo Methods in Statistical Physics 2.Application of Monte Carlo Method in Statistical Physics
