分子模拟教程

第二章数值积分和MonteCarlo方法第一节数值积分baSfxdx令10,,kknhxxxaxb，则110(),''()kknkkkkkkkkxxSfxdxfxfxxfffxffxxfakxbx零阶近似hfxfk1010nkknkkhfhhfhS一阶近似21hhffxxfxfkkkk∵1212212122kkxxkkkkkkhxxxxxdxxx∴102121nkkkkhhffhfS102121nkkkhffh从直观看，用112kkff近似fx比只用kf或1kf好。这方法也称Trapezoid方法。这样的数值积分方法的优点：简单直观，误差可以控制缺点：“平均主义”，在0xf的区域，kfxx对S贡献很小，但消耗同等的机时。在多自由度系统这弱点尤为特出。问题：直观地看，零级近似和一级近似的差别在哪？习题：编程序数值计算高斯积分。第二节MonteCarlo方法如何用随机方法求积分？例如，可用‘抛石子’方法。但这方法不比简单的数值积分有效。1．简单抽样的MonteCarlo方法均匀地随机地选取[ba,]中kMx个点，显然，11()1/MkkSfxMM当M足够大，当然可以得到足够好的积分值。问题：为什么误差是1/M？答：不妨把这看成一个M次测量的实验，假设每次测量都是独立的，由涨落理论，误差应为1/M。比较误差：MonteCarlo方法1/SM数值积分Trapezoid方法21/(1/)SMhM对单自由度而言，数值积分方法要有效得多。对多自由度，例如d自由度，MonteCarlo方法1/SM！！数值积分Trapezoid方法2/1/(1/)ddSMhM当d非常大，数值积分方法根本没法和MonteCarlo方法比较。我们当然可以再改进数值积分方法的精度，但这种改进的量级没法和d的大小比拟。在多体系统的数值模拟中，d通常至少是410！MonteCarlo方法真的是完美了吗？当然不是。‘平均主义’的弱点其实还没改进下面我们引入所谓的重要抽样MonteCarlo方法当引入重要抽样方法后，每次抽样的样品可能不独立如何取得独立的抽样，是MonteCarlo方法的重点所在！2．重要抽样的MonteCarlo方法如果被积函数f(x)是均匀的函数，则简单抽样方法已经可以得到相当准确的积分值。如果被积函数f(x)不是均匀的函数----这在高维积分十分常见，则必须引入重要抽样。我们希望在对积分贡献大的区域多取样品，在贡献小的区域少取样品。设积分为()()baSfxfxWxdx其中f(x)是均匀函数，W(x)是不均匀函数，而且()0Wx，可以归一化，给予概率分布的含义。假设我们可以按照分布W(x)得到{}lMx个点，则11()1/MllSfxMM注意：现在W(x)不出现在求和式子中，而是体现在{}lx的分布里。证明：如第一节方法，把[ba,]分为n个小区间，,lkkxxxh落在区间的数目约为kWxhM对,lkkxxxh1011lkMnlkklkfxfxfxfxWxhMMM显然1limlimMnlkkMnlkfxfxWxhfxSM关键：如何按分布Wx，产生M个{}kx（即上面的{}lx）问题：可以用产生随机数的方法产生Wx吗？答：多自由度，有相互作用时不可能3．Markov过程产生{}kx的方法：构造一个Markov过程，即给出一个动力学规则，由tx随机地产生1tx。那么，给一个0x，可产生0001,,,MMMxxxx如果随机过程满足一定条件，则当0M足够大时，即达到“平衡态”时，001,MMMxx按Wx分布。两个条件：各态历经从概率上说，在有限时间内tx可走遍[ba,]细致平衡Markov过程由从tx到1tx的转移概率定义。设Txx为过程从x转移（或跃迁）到x的概率，Txx为从x到x的概率。则细致平衡条件为TxxWxTxxWx记忆：从x转移到x的概率正比与Wx证明：设,Wxt为x的“非平衡态”分布,,,xxdWxtTxxWxtTxxWxtdt,,,0tWxtWxdWxtdt,,xxTxxWxTxxWx细致平衡条件显然保证Wx满足上述方程，所以，,WxWx不过，细致平衡条件是充分条件。第三节Metropolis算法和Heat－bath算法Markov过程的全部信息包含在转移概率Txx中细致平衡条件是平衡态Wx的要求是否各态历经常常可以直观判断Txx必须给出从任意x到任意x’的概率，但这并不一定要求从x到任意x’的概率都是零。各态历经只要求在有限时间内，即Txx的有限次作用，能到达x’。1．Metropolis算法Metropolis(1915-1999)ThePaperwascited7500timesfrom1988to2003令(')/()1(1,(')/())mTxxWxWxWxWxWxWxMinWxWx设WxWx(')/()(')1()mmTxxWxWxWxTxxWx设WxWx1(')()/(')()mmTxxWxTxxWxWxWx即mTxx满足细致平衡条件。但是，注意到Txx是转移概率，所以应有归一化条件'1xTxx上面的mTxx还不满足归一化条件。所以，完整的Metropolis算法的转移概率是()mTxxSxxTxx其中()Sxx是从x选中x’的概率，而且()(')SxxSxx。这样的转移概率仍然满足细致平衡条件。当然，归一化条件也能保证。例如，可选1．()Sxx常数，即均匀选取x’2．'[,]0'[,]xxxSxxxxx常数归纳起来，Metropolis算法包括如下步骤：1）设txx，按()Sxx选取尝试x’，2）以概率mTxx取1'txx以概率1mTxx取1txx（！！）第二步要记住，初学者容易误解，以为一旦x’不被采纳，重新选取尝试x’。不难理解，无论是1或2方案的()Sxx，Metropolis算法都满足各态历经条件，因为只要(')0Wx，在有限时间取到x’的概率不为零。Metropolis算法非常普遍。一个重要的例子，在物理学中，常取()Wx为正则分布，()/()HxkTWxe其中H(x)代表能量，T是温度，k是Boltzmann常数。所以，跃迁概率为//1(1,)HxHxHxHxkTmkTHxHxHxHxeTxxMine如何在计算机上实现？①设已到达txx点，按均匀分布随机地取,,xxx为一‘恰当’的数例如，取1,0r为计算机上常用的“随机”数，则21xxr②以概率mTxx取1'txx以概率1mTxx取1txx例如，计算xHxH，取随机数1,0r如果/HxHxkTer，取1'txx，否则1txx/HxHxkTe③如何取0MM和？这与具体系统有关。通常通过尝试改变0M以判断0M是否够大，即随机过程是否已到达“平衡”态。一般010MM。不过，对单自由度问题，0M不重要。练习：①221xxH计算2x，和解析结果比较②424121xxxH，计算2x讨论的作用对没有解析解的系统，如何判断数值结果的可靠性和正确性？有兴趣的同学还可以和数值积分比较，会发现MonteCarlo方法在单自由度情形没有优势2．Heat-bath算法Heat-bath算法没有Metropolis算法那么普遍但也是多体系统研究中的典型方法之一Heat-bath算法直接取转移概率为(')hTxxWx这算法自然是各态历经并且满足细致平衡条件的。Heat-bath算法的特点是hTxx与x无关在某些情形会显示出优越性当(')Wx较复杂，在计算机上不一定能够找到简单的实现事实上，对单自由度情形，已经是直接产生平衡态分布()Wx例如，如果x只取x1和x2,可取(1)'1(2)'2hWxxxTxxWxxx计算机上的实现：取随机数1,0r如果(1)Wxr，取11txx，否则12txx这一算法可简单推广到x取多个分立值情形。当x是连续变量时，不一定能实现Heat-bath算法。但我们可以结合Metropolis算法和Heat-bath算法。例如，取1()hTxxSxxTxx1111()/(()())'()/(()())'0hWxWxWxxxTxxWxWxWxxxotherwise这算法Metropolis的特征多些。归纳起来，Metrpolis算法或Heat-bath算法可以相当普遍地实现重要抽样的MonteCarlo模拟方法。但是，对多体系统，抽样得到的样品可能不独立，即有关联。设关联时间为,即每个样品才有一个是独立的。那么，误差为(/)SOM当很大时，MonteCarlo方法遇到极大困难。对物理学家而言，这是MonteCarlo算法需要解决的重要问题。第四节Langevin方程问题：系统的Hamiltonian量为S，平衡态分布为)exp(S，（这里温度已吸收到S）。假设系统0t时处于一初始状态，系统如何演化至平衡态？单自由度实数:xxSSLangevin方程txxSdttdx:02ttttt高斯随机数对固定tP2e2201edZedZt如果没有随机力，平衡态为0dxtSxdtx，即能量取极小值。如果存在随机力，体系会被推离能量极小，处于某种能量较高的平衡态。由于随机力的存在，Langevin方程存在复杂性，因为我们必须考虑对随机力带来的奇异性。为了简单起见，我们对时间分立化在数值模拟中应用较直观ttttt2，tttttt01Z21=21lnZ∴4tLangevin方程()Sxtxtxttxttttx令2t∴()(')ttttttttxtxStx2方程的解tx是随机变量，在数值模拟中给定初始值txx,0还不确定，与随机力有关。也就是说，在t时刻，x遵从一个分布;Pxt。物理量x的平均值;xtdxPxtx问题：xt的含义？答：必须对t之前的所有随机力做平均。Langevin方程其实也是一种MonteCarlo方法设')(,)(xttxxtx]4/))('(exp[]2/)(exp['22ttxxSxxtxxP]4/)')'('(exp[]2/)(exp['22ttxxSxxtxxP注意到x和x’之差是无穷小，)]()'((exp[])()