博士生讲座2015解析.ppt

蒙特卡洛方法经历了几十年的发展，基本算法有很多，其中最为著名的算法是Metropolis等人提出的通过马尔科夫链（Markov Chain）的算法进行系统抽样。使用MC方法进行分子模拟计算时，马尔科夫链从某一初始状态出发，使用随机数发生器产生一系列的随机的分子构型，以粒子移动的计算过程为例： 式中Δmax为给定的最大位移，ranf()为[0,1]上符合均匀分布 的一个随机数。如果采用构型偏倚MC方法进行模拟，则生成的随机数不是在[0,1]区间均匀产生的。 但是这样产生的初始构型并不一定满足Boltzman分布，所以为了保证该初始构型能被接受，必须要满足细致平衡：在某确定体系中，从空间某一点o向另一点n移动的平均可接受数目K(o→n)，应等于空间中点n向点o移动的平均可接受数目K(n→o)。 （1.15） 上三式中，N(o)和N(n)分别表示体系中点o和点n存在粒子的概率；α(o→n)和α(n→o)分别表示在体系中，粒子从点o向点n移动以及粒子从点n向点o移动的概率；和分别表示从粒子点o移动到点n以及粒子从点n移动到点o的接受概率。在Metropolis抽样方法中，α(o→n)=α(n→o)，而在构型偏倚抽样方法中，α(o→n)≠α(n→o)。当新的构型r在体系中平衡时，体系中的配分函数表示为Z，表达式如下： 其中U(rN)为新构型下体系的能量，构型r存在的概率密度为： 为满足能量最低原理，必须满足下面的条件：如果U(n)U(o)，则需要将exp[-β(U(n)-U(o)]与一个[0, 1]区间的随机数进行比较，如果exp[-β(U(n)-U(o)]小于该随机数，就不接受本次移动，仍将原构型作为马尔科夫链的初始状态，进行新的MC模拟；如果exp[-β(U(n)-U(o)]大于该随机数，则接受本次移动。当U(n)≤U(o)时，则本次移动被接受，新的构型作为马尔科夫链的一个新的初始状态进行下一步的迭代。 分子动力学模拟（MD) 分子动力学模拟是基于牛顿力学基础上一种分子模拟方法，通过求解体系所有粒子的运动方程，得到各个时刻每个粒子的坐标与速度，不同时刻的坐标与动量数据描述了体系在不同时刻的微观状态。体系的宏观性质通过对微观状态的时间平均获得。 最早的MD模拟是在1957年，由Alder和Wainwright 实现的硬球势模型，早期分子动力学模拟在时间尺度和空间尺度都非常小。上个世纪80年代以来，计算机技术的发展和统计系综技术的完善，分子动力学模拟逐渐走向了实用阶段。 为讨论简单，考虑体系由N个全同粒子组成，粒子放置在体积为V、空间维数为d 的容器内，体系的温度为T。i粒子的坐标 和动量 满足经典力学的Hamilton（哈密顿）正则运动方程 （7） 式中，H为体系的Hamilton函数，即粒子体系的动能和势能之和。 如若考虑粒子间的相互作用为两体势u(rij),Hamilton函数H则为 其中m为粒子的质量，u(rij)是分别处于ri和rj的两粒子间的相互作用能，rij=|ri-rj|，正则方程具体化为 （8） （9） 上式中的 为梯度算符，合并上两式，得 式（10）左边是第i个粒子的加速度，方程的右边是第i个粒子受到其他粒子的作用力之和除以粒子的质量，表示该方程就是第i个粒子所满足的Newton运动方程。 （10） 在势函数、初始条件和边界条件给定后，原则上就可以利用计算机求解粒子运动方程组了。 MD模拟中的关键步骤就是解粒子运动方程，获取体系在不同时刻的微观状态。计算机解微分方程的基本方法是将时间和空间离散化，用差商代替微商，将微分方程简化为差分方程。MD模拟中的力的计算是最耗时的步骤。已发展了多种有限差分方法。Verlet算法，Leap-frog算法和Velocity-Verlet算法 分子模拟法的计算量要比量化计算小得多，而且能应用于几百乃至上万个原子的系统。而且在适当的范围内，分子模拟方法的计算精度与量子化学计算相差无几。对大分子复杂体系而言，分子力学方法是一套行之有效的方法。 分子模拟需要知道粒子之间的相互作用力