《化学统计热力学》课件Statistical Thermodynamics7

1、Chapter 7 Monte Carlo及相关方法7.1 Monte Carlo模拟7.2 分子动力学模拟7/11/2022第七章 Monte Carlo方法27.1 Monte Carlo模拟派什么用场；用计算机来模拟统计的随机抽样(Sampling)过程有非常广的应用，只要能把问题转化为随机抽样，即可应用应用于事件的几率已知，但结果很难预知的过程名字的由来Monaco（摩纳哥）公国的一个城市，以赌博闻名。从大约1944年开始发展和完善7/11/2022第七章 Monte Carlo方法37.1 Monte Carlo模拟7/11/2022第七章 Monte Carlo方法4Monte C

2、arlo的要素Probability distribution functions (pdfs)研究的对象必须可以转化为一些几率分布函数来描述随机数(Random Number)产生均匀分布的等间隔的随机数Sampling rules （抽样规则）如何取样（物理上的测量如何进行）7/11/2022第七章 Monte Carlo方法5Monte Carlo的要素, continued.如何评价(scoring)感兴趣的量把每次取样综合起来得到感兴趣的量的值（物理量的平均值）误差估计与随机取样的样本数目或是其它量的关系尽量减少变量数目以节省Monte Carlo模拟的时间并行化的考虑蒙特卡罗方法的

3、收敛性，误差 蒙特卡罗方法作为一种计算方法，其收敛性与误差是普遍关心的一个重要问题。收敛性误差减小方差的各种技巧 效率收敛性 蒙特卡罗方法是由随机变量X的简单子样X1，X2，XN的算术平均值： 作为所求解的近似值。由大数定律可知，如X1，X2，XN独立同分布，且具有有限期望值（E(X)），则 即随机变量X的简单子样的算术平均值 ，当子样数N充分大时，以概率1收敛于它的期望值E(X)。误差 蒙特卡罗方法的近似值与真值的误差问题，概率论的中心极限定理给出了答案。该定理指出，如果随机变量序列X1，X2，XN独立同分布，且具有有限非零的方差2 ，即 f(X)是X的分布密度函数。则 当N充分大时，有如下

4、的近似式 其中称为置信度，1-称为置信水平。 这表明，不等式 近似地以概率 1成立，且误差收敛速度的阶为O(N-1/2) 。 通常，蒙特卡罗方法的误差定义为上式中与置信度是一一对应的，根据问题的要求确定出置信水平后，查标准正态分布表，就可以确定出。几个常用的与的数值：关于蒙特卡罗方法的误差需说明两点：第一，蒙特卡罗方法的误差为概率误差，这与其他数值计算方法是有区别的。第二，误差中的均方差是未知的，必须使用其估计值来代替，在计算所求量的同时，可计算出 0.50.050.003 0.67451.963减小方差的各种技巧 显然，当给定置信度后，误差由和N决定。要减小，或者是增大N，或者是减小方差2。

5、在固定的情况下，要把精度提高一个数量级，试验次数N需增加两个数量级。因此，单纯增大N不是一个有效的办法。另一方面，如能减小估计的均方差，比如降低一半，那误差就减小一半，这相当于N增大四倍的效果。因此降低方差的各种技巧，引起了人们的普遍注意。后面课程将会介绍一些降低方差的技巧。 效率 一般来说，降低方差的技巧，往往会使观察一个子样的时间增加。在固定时间内，使观察的样本数减少。所以，一种方法的优劣，需要由方差和观察一个子样的费用（使用计算机的时间）两者来衡量。这就是蒙特卡罗方法中效率的概念。它定义为其中c是观察一个子样的平均费用。显然 越小，方法越有效。 7/11/2022第七章 Monte Ca

6、rlo方法13一般程序产生随机数转变为随机变量（抽样）利用变量来模拟产生结果Analyze Output是否继续7/11/2022第七章 Monte Carlo方法141.随机数的产生真正的随机数是很难产生的伪随机数(Pseudo Random Number)产生的方法丢骰子，掷硬币摇奖做好的表格计算机递推方法7/11/2022第七章 Monte Carlo方法15对随机数发生器(Random Number Generator, RNG)的要求产生0,1间均匀分布的随即序列各随机数间是独立的数目巨大，很少重复可重复得到有效，并且容易移植7/11/2022第七章 Monte Carlo方法16R

7、NG-乘同余法Zi = mod(aZi-1+b,m) m: 通常是计算机可以表示的最大正整数，如(231 1)a 正的奇数，如 75 = 16807b 步长，整数Zo (Starting) = Zi-1 -初值，种子7/11/2022第七章 Monte Carlo方法17RNG-乘同余法随机数i注意，取a和m时须小心，避免循环7/11/2022第七章 Monte Carlo方法18乘同余法的改进：混合伪随机数序列乘同余法的缺点：伪随机数之间的相关性比较强改进混合伪随机数首先产生一组整数随机数列，然后从其中取出一个整数，根据这个整数再产生新的随机数列结果相当于将一个伪随机数列按另一个伪随机数列的

8、大小重排，因而减弱了所用的随机数序列的相关性。缺点：增加计算机时7/11/2022第七章 Monte Carlo方法192.把随机数变为变量随机数与随机变量要一一对应逆变换的方法产生一个随机数由逆函数X=F-1()求得变量X011020300.77/11/2022第七章 Monte Carlo方法20分布函数种类均匀分布指数分布高斯分布三角形分布7/11/2022第七章 Monte Carlo方法21均匀分布(U-L)+L : 随机数L: 变量的下限；U: 变量的上限例：在公共汽车站台某一路车到达的时间7/11/2022第七章 Monte Carlo方法22三角形分布7/11/2022第七章

9、Monte Carlo方法23指数型分布如Boltzmann分布7/11/2022第七章 Monte Carlo方法24高斯型（正则）分布如考试成绩7/11/2022第七章 Monte Carlo方法25例子：对函数的积分定域积分：产生随机数， 0,1抽样：得到随机变量f()，根据模拟：得到感兴趣的量7/11/2022第七章 Monte Carlo方法26例子：对函数的积分7/11/2022第七章 Monte Carlo方法273. 抽样方法：理论上U可取到无穷大，实际很大的值对配分函数的贡献值为0分布函数重要性实际体系配分函数的构型积分e-U/kTU7/11/2022第七章 Monte Ca

10、rlo方法28重点抽样法G(x)满足G(1)=1, G(0)=0，将G(x)看作是区间0,1上的分布函数，则g(x)为几率密度函数有目的的抽样减小方差，加速收敛7/11/2022第七章 Monte Carlo方法29重点抽样法，continued积分值方差如7/11/2022第七章 Monte Carlo方法30重点抽样法，continuedBoltzmann分布的重点抽样7/11/2022第七章 Monte Carlo方法31Boltzmann分布的重点抽样，continued现在，分布函数变成了7/11/2022第七章 Monte Carlo方法32Boltzmann分布的Metropol

11、is抽样方法思想：是抽样尽可能按实际（也就是平衡态）的分布进行做法：构建一个Markov链，每改变一个构型，即改变体系中分子的坐标pi,qi，改变一次，即为走一步(A Step)每走一步，得到相互作用能(势能）U，即判断是否合理，保留合理的走法求物理量的期望值7/11/2022第七章 Monte Carlo方法33Metropolis抽样方法产生随机数0,1改变每个分子的坐标步长计算得到体系的势能Unew及势能变化Unew半经验的方法计算势能，Lenard-Jones势函数精确从头算计算7/11/2022第七章 Monte Carlo方法34Metropolis抽样方法, continue另外

12、再产生随机数rand(0,1)判断true，保留改变，继续对新坐标改变false，重新改变构型完成足够数量的循环次数后，得到体系物理量（热力学函数）的平均值7/11/2022第七章 Monte Carlo方法35Metropolis抽样方法, continue步长rmax的注意事项不能太大，否则很多次改变不合理不能太小，否则效率差，需走很多步才能探索完整个构型空间7/11/2022第七章 Monte Carlo方法364. 分子的Monte Carlo模拟比原子体系复杂改变分子质心坐标的同时还要考虑分子的旋转不成功的step更多刚性分子只须考虑分子的整体旋转绕x轴选旋转角的坐标变换7/11/2

13、022第七章 Monte Carlo方法37分子的Monte Carlo模拟，continue任意旋转，Euler angles7/11/2022第七章 Monte Carlo方法38分子的Monte Carlo模拟，continue任意旋转的坐标变换实际上对cos的做法是：7/11/2022第七章 Monte Carlo方法39分子的Monte Carlo模拟，continue使用Euler角的缺点上面的矩阵含有6个三角函数，计算机计算耗时解决办法：变量替换7/11/2022第七章 Monte Carlo方法40分子的Monte Carlo模拟，continue坐标变换矩阵即变为：因此需要产

14、生四个随机数实现对q0,q1,q2,q3进行改变：产生-1,1间的随机数1, 2，使S1= 12+221；产生-1,1间的随机数3, 4，使S2= 32+421；变换的基失(q0,q1,q2,q3)为：7/11/2022第七章 Monte Carlo方法41柔性分子或聚合物高分子的模拟分子的内部自由度太多，非常困难固定某些内部自由度，如键长和键角，而放开两面角高分子的模拟请参考相关专著，如：杨玉良, 张红东 著：高分子科学中的 Monte Carlo 方法，复旦大学出版社，19937/11/2022第七章 Monte Carlo方法425. “Biased” Monte Carlo Metho

15、dsExplore the most important part of the phase spaceSingle solute molecule surrounded by a large number of solvent moleculesPreferential sampling接近溶质分子的溶剂分子移动比离得远的快7/11/2022第七章 Monte Carlo方法43Preferential sampling选取一个区间(cutoff region)，在这个区间外的分子移动慢而区间内的分子移动快选取一个几率p0，对一个分子，如果在cutoff region内则移动反之，产生一个0

16、-1间的随机数p1p1 p0，不动Cutoff region的选取也可以不固定7/11/2022第七章 Monte Carlo方法44Preferential sampling, continue需要注意的问题移出cutoff区间内的分子比移进来的多Force biased Monte Carlo原子或分子按照它所受的力的方向移动Calculate forcesA: 常数fi: 力riG: 按Gaussian分布，且平均值为0，variance为：2A7/11/2022第七章 Monte Carlo方法456. 不同系综的Monte Carlo方法最简单的就是正则系综即NVT系综NPT系综：改

17、变原子或分子坐标的同时改变体积体积改变后，体系的势能也要重新计算，对于简单的相互作用势能，下面的技巧：scaled coordinates是很有用的（Lold是体系所处立方体的边长）:7/11/2022第七章 Monte Carlo方法46NPT系综上式可以重新写为：体积改变后体系的新势能为：7/11/2022第七章 Monte Carlo方法47NPT系综, continue则势能变化：对于NPT系综，判断一个step要不要继续是根据H：7/11/2022第七章 Monte Carlo方法48巨正则系综(VT系综）化学势可写作z：活度，常数de Broglie 波长7/11/2022第七章

18、Monte Carlo方法49巨正则系综, continue需要用到随机数作改变的量(取样)：粒子的运动，用Metropolis方法取样粒子的湮灭(destroy)粒子的生成(create)7/11/2022第七章 Monte Carlo方法50巨正则系综, continue一些注意事项粒子产生的几率和粒子湮灭的几率要一样对于一些比较密集的体系，产生和湮灭成功的几率比较小7/11/2022第七章 Monte Carlo方法517. 计算化学势按照化学势的定义N很大7/11/2022第七章 Monte Carlo方法52化学势，continue对于实际的气体或液体，气体或液体的构型积分：单个分子

19、（理想气体）配分函数7/11/2022第七章 Monte Carlo方法53化学势，continue则化学势，令，7/11/2022第七章 Monte Carlo方法54化学势，continue化学势，7/11/2022第七章 Monte Carlo方法55: Configurational bias Monte Carlo Method参考：Siepmann, J. I. A Method for the Direct Calculation of Chemical Potentials for Dense Chain Systems. Mol. Phys. 1990, 70(6), 1145-1158.Siepmann, J. I.; Frenkel, D. Configur