中国矿业大学分子模拟课件第十二章.pdf
《中国矿业大学分子模拟课件第十二章.pdf》由会员分享,可在线阅读,更多相关《中国矿业大学分子模拟课件第十二章.pdf(43页珍藏版)》请在三一文库上搜索。
1、分 子 模 拟 牛继南 2011.4 第十二章 蒙特卡罗计算方法 12.1 基本原理 广义定义:指凡是在计算中引用随机数的 计算法,英文名称Monte Carlo,简称MC。 以计算半径为单位长度的圆的1/4面积为 例说明蒙特卡罗方法基本原理。 选取一组随机点P(x,y),它们的范围为 0x1,0y1。 如果此点满足 则被接受,否则被排除。 如果随机点数为N,被接受的点数为n,则 我们利用随机点求得的1/4圆面积为: P点落入1/4圆面积的概率乘以11正方形面 积。 1/4圆的真实面积为r2/4=0.785398,而不 同随机点数得到的面积为: 以上的问题相当于用蒙特卡罗法求积分: 可以将这
2、种方法推广求解多重积分,一般 用于普通方法难以求解的积分,如对于: 计算时,在变量xi可能的区间内任意取n组 随机点 ,i=1,n. 相当于计 算: 这种方法的误差为n-1/2,故选取的点越多, 计算的积分值越精确。 对于积分函数: 可以改写为: 其中p(x)为配分函数,设在范围x1,x2内按 照p(x)分布选取n个随机数x,则依照统计学 积分式可写为: 式中表示对n个随机数的平均。 对于正则系综中任意坐标,函数F的平均值 为: 式中P(x)为: 为系统位于特殊结构 的概率。 定义结构积分为: 则F的平均值为: 理论上若任意选取大量的结构点,则可按 MC法求解积分,一般步骤为: 1. 随机产生
3、N个原子坐标; 2. 计算此构型的势能和函数F; 3. 计算玻尔兹曼因子 4. 累加玻尔兹曼因子,并将函数与玻尔兹曼因子相 乘后累加,回复步骤1; 5. 当计算Nt循环后,求出F的平均值为: 理论上可以用蒙特卡罗方法计算统计的平 均值,但若计算大体系必须选取大量的样 品点,而这些样品点所对应的结构多数因 为能量高,因此不容易出现。此外,依靠 此方法计算积分时,必须先计算结构积分 ,但实际上此积分不易收敛。 基于这些原因,事实上直接用MC法并不实 际,通常采用Metropolis的重要取样法。 12.2 重要取样法 以NVT系综的某一物理量A的平均值为例, 统计力学算式为: 式中NVT为概率,
4、为对应的结构。从统计 角度讲,有些结构出现的概率大(trial), 有些结构出现的概率小,而重要取样法计 算平均时设法选取重要的 ,使积分仅为 这些重要贡献的求和: n为重要取样点。 这种方法是Metropolis,Rosenbluth,Teller 等在1953年首先提出,被称为Metropolis法 ,是现在MC计算的基础。 此法的重点是建立系统的马科夫链(markov chain),使最终满足重要的取样点的要求. 马科夫链符合两个条件: 1. 设系统个状态所成的有限集合为 则每次选择的随机数取样点产生的状态 必为此 集合的一个状态。 2. 每次所选的取样点仅与前一个取样点有关,而与 由何
5、处开始取样无关。 以计算机运行为例: 若当天正常运转,则次日正常运转的概率 为60%;若当天停机,则次日也停机的概率 为70%。 考虑此问题时分为正常(状态1)和停机( 状态2)两种状态,由一种状态转换成另一 种状态的概率用ij来表示,则: 11表示第一天正常第二天也正常,为0.6 12表示第一天正常第二天停机,为0.4 21表示第一天停机第二天正常,为0.3 22表示第一天停机第二天也停机,为0.7 明显有: 可以建立转换矩阵 设1为系统处于状态1的概率,2为系统处 于状态2的概率,则系统的状态概率可用概 率向量表示: 设 为系统起始状态的概率向量, 为系 统经过第一次转换后的概率向量,即:
6、 再经过一次转换后: 以此类推,即得到马科夫链。 如果第一天计算机正常运转和停机的概率 相等,即 =(0.5,0.5),则第二天的概率 向量为: 即第二天正常概 率为0.45,停机 概率为0.55 第三天的概率向量为: 经过多天的计算得到极限值: 若概率达到极限,则再乘转换矩阵必得到 相同的概率向量,如: 这时表明系统的概率分配值不再变化,系 统达到平衡。 推广到N中可能的系统,则 为N维向量, 转换矩阵 为NN矩阵,矩阵中每行元素 的和为1,即 满足极限概率值的条件为: 为了得到这个方程的解,通常采用微观可 逆条件来替换,即: 如果对m的所有态求和,则 = n m nmnnm m nmn m
7、 m 即为上面的式子。 第一个提出解这个方程的为Metropolis等, 他们引入对称矩阵,即mn=nm,表示选择 mn和选择nm的概率相同,设 分 别为状态m和状态n被接受的极限概率, 若已达到极限状态,且 ,则m转换为n 的概率等于选择mn的概率,即 若 则m转换为n的概率等于选择 mn的概率乘以n状态与m状态被接受的概 率比。 同时还需满足同一状态的转换概率: 由于多粒子系统有不同的状态,故转换矩 阵十分庞大,但极限状态可以由统计力学 得到。对任一结构 ,其极限概率为 因此 式中, 为二状态n与m的势能差。 当n状态的势能低于m状态的势能,则此结 构的转换可以被接受,若n状态的势能高于
8、m状态的势能,则根据前面的式子来接受此 移动。 实际当中可将玻尔兹曼因子与一个在0-1之 间的随机数比较,如果大于次数则被接受 ,如果小于次数则排除此转换。 按照这个图,随机数在左边靠近坐标为0的 时候,配分函数值较大,容易出现在函数 下方,所以容易被接受(此时虽然UnUm ,但能差很小);如果随机数在右边,配 分函数值接近于0,因此随机点易出现在函 数上方,被排除(此时UnUm,且能差很 大)。 因此通过这样的转换可排除不可能的移动 ,达到重要取样的作用。 转换过程中转换矩阵 随系统的不同而不 同。 12.3 Metropolis具体计算步骤 以二维空间为例,系统中有数 个粒子,处于系统m。
9、产生另一 状态时随机选择系统中的一个粒 子i,将其由原来的位置 以相同 的概率移动到一定范围R内的任一位置 。 范围R的中心为 ,边长为 。这样可以 产生NR个新位置,即NR个新的状态,因此 选择mn的概率为1/NR,可定义 为: 下一个步骤为计算移动的玻尔兹 曼因子。设选择的粒子i由原来的 移动至新的位置 ,设系统中 的粒子作用仅为LJ势,则: 即只需计算能量变化量。若 ,将玻 尔兹曼因子和随机数 比较,若大于 次数则被接受,反之则被拒绝。若接受以 新的位置为i原子位置,否则回到原位。 对含有N个粒子的NVT系统,Metropolis的 计算步骤如下: 1.任一挑选一个粒子,产生任意方向的位
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 中国矿业 大学 分子 模拟 课件 第十二
链接地址:https://www.31doc.com/p-3698160.html