分子模拟教程径向分布函数

分子模拟教程径向分布函数第一页，共七十三页，编辑于2023年，星期五掌握分子模拟方法的必备知识：编程技能(FortranorC/C++)

统计物理学（统计力学）：

统计物理学基础；系综原理；

非平衡统计力学基础；涨落理论分子热力学

：

分子间相互作用理论；分布函数理论气体分子运动论其它第二页，共七十三页，编辑于2023年，星期五分子模拟的目的:为什么要进行分子模拟？将分子聚集体的性质与如下方面相联系：分子的微观相互作用分子聚集体的结构分子的动力学过程分子模拟对实验进行补充，使我们能够：预测现有或新材料的性质在分子水平研究宏观现象获得实验无法或难以发现的东西第三页，共七十三页，编辑于2023年，星期五什么是计算机分子模拟方法？分子模拟的定义：

统计力学基本原理出发，将一定数量的分子输入计算机内进行分子微观结构的测定和宏观性质的计算。按照获得微观态的方法不同，分子模拟分为：蒙特卡罗方法(MonteCarlo,MC)分子动力学方法(MolecularDynamics,MD)(3)混合方法(hybridmethod，HM)第四页，共七十三页，编辑于2023年，星期五计算机分子模拟的发展历史：1.蒙特卡罗方法（MC）

1953 Metropolis,Ulam,RosenbluthandTell，

LosAlamosNationalLab MonteCarlosimulationofhardsphere.2.分子动力学方法（MD）

1957 AlderandWainwrigth， LivermoreLab

Moleculardynamicssimulationofhardspheres.第五页，共七十三页，编辑于2023年，星期五微观与宏观分子模拟在微观尺度与实验室的宏观世界之间起着桥梁的作用：给定分子间的相互作用

“准确”预测研究体系的性质第六页，共七十三页，编辑于2023年，星期五MC与MD的区别：MC:

构型平均，不包含动力学部分；

利用概率行走产生微观态。

MD:时间平均，产生动力学性质；

利用运动轨线随时间的变化来产生一系列微观态。第七页，共七十三页，编辑于2023年，星期五计算机分子模拟的发展历史（续）：从上个世纪九十年代初期以来，计算机模拟技术得到了飞速发展，主要基于三个方面的发展：

分子力场的发展（基石）

（Amber，OPLS、Compass）

原子间的键长、键角、分子间的内聚能等

模拟算法（途径）

计算机硬件（工具）HPCx第八页，共七十三页，编辑于2023年，星期五计算机分子模拟的特点：原子水平的模拟计算机实验检验理论、筛选实验科学研究中的第三种方法第九页，共七十三页，编辑于2023年，星期五分子模拟中涉及的几个基本概念：模拟计算盒子或模拟胞腔Simulationbox(cell)装有一定数目流体分子的研究对象，它是我们要研究的宏观体系的缩微模型。立方形胞腔第十页，共七十三页，编辑于2023年，星期五周期性边界条件(Periodicboundarycondition,PBC)在小体系中，边界效应总是很显著。在包含1000个原子的简单立方晶体中－488个原子处于边界上。在包含1000000个原子的简单立方晶体中－仍然有6%的原子在边界上。在模拟中，考虑具有真实边界的对象，不切合实际：增强了有限尺寸效应人为造成的边界会影响流体的性质第十一页，共七十三页，编辑于2023年，星期五当某个粒子运动出模拟盒子的某一边界时，另外一个影像粒子从另一对立边界进入到此盒子中。周期性边界条件(Periodicboundarycondition,PBC)本体体系的近似：中心盒子在X，Y和Z方向无限扩展；消除人为形成边界的表面效应；保证中心盒子中的粒子数恒定。只需要跟踪中心盒子中各粒子的运动。第十二页，共七十三页，编辑于2023年，星期五周期性边界条件的算法：yxLo采用数学函数：FLOOR(r/L):返回不超过r/L的最大整数FLOOR(4.8)hasthevalue4.FLOOR(-5.6)hasthevalue-6.第十三页，共七十三页，编辑于2023年，星期五xL/2o-L/2采用数学函数：r/L>0,ANINT(r/L)=AINT(r/L+0.5)r/L0,ANINT(r/L)=AINT(r/L-0.5)周期性边界条件的算法：y第十四页，共七十三页，编辑于2023年，星期五最小影像转化原理(Minimumimageconvention)定义：中心元胞中的一个粒子只与此元胞中的其它N－1个粒子，或它们的最近邻影像发生相互作用。适用条件：粒子间相互作用势能的截断距离必须不大于模拟中心元胞长度的一半。此两粒子与中心粒子的距离相等，但是：黑色球发生作用

绿色球不发生作用此两粒子是与中心原子相互作用的最近邻影像第十五页，共七十三页，编辑于2023年，星期五最小影像转化原理的算法：采用数学函数：r/L>0,ANINT(r/L)=AINT(r/L+0.5)r/L0,ANINT(r/L)=AINT(r/L-0.5)第十六页，共七十三页，编辑于2023年，星期五截断势能(TruncatingthePotential)本体体系采用周期性边界条件描述：不可能将所有粒子与它们影像粒子间的相互作用全都计算。必须在不大于中心盒子长度的一半处进行截断，以便与最小影像转化原理一致。粒子间的相互作用主要来自于截断范围内，而范围外的贡献很小，可忽略不计。截断范围内的相互作用第十七页，共七十三页，编辑于2023年，星期五截断势能函数的形式：简单截断势能函数(TruncatedPotential)：缺点：rc:截断距离或半径势能在截断处不连续,当一对分子穿越边界时，总能量不守恒。分子间力在截断处为无穷大，MD运动过程不稳定。忽略截断半径之外的所有作用第十八页，共七十三页，编辑于2023年，星期五位移截断势能函数(ShiftedandTruncatedPotential)：缺点：分子间力仍然在截断处不连续。优点：势能在截断处连续，但不影响分子间力的大小分子间力在截断处不为无穷大截断势能函数的形式：常用于MC和MD模拟中第十九页，共七十三页，编辑于2023年，星期五位移－力截断势能函数(Shifted-ForcePotential)：常用于MD模拟中优点：势能和分子间力均在截断处连续截断势能函数的形式：第二十页，共七十三页，编辑于2023年，星期五截断势能函数的对比：位移－力截断势能简单截断势能函数第二十一页，共七十三页，编辑于2023年，星期五一、MonteCarlo模拟方法基础：亦称统计模拟或随机抽样方法，statisticalsimulationmethod

利用随机数进行数值模拟的方法MonteCarlo名字的由来：是由Metropolis在二次世界大战期间提出的：Manhattan计划，研究与原子弹有关的中子输运过程；NicholasMetropolis(1915-1999)Monte-Carlo,Monaco投硬币，掷骰子第二十二页，共七十三页，编辑于2023年，星期五MonteCarlo方法计算Pi值第二十三页，共七十三页，编辑于2023年，星期五随机数的定义和特性什么是随机数？单个的数字不是随机数;是指一个数列，其中的每一个体称为随机数，其值与数列中的其它数无关；在一个均匀分布的随机数中，每一个体出现的概率是均等的；例如：在[0,1]区间上均匀分布的随机数序列中，0.00001与0.5出现的机会均等第二十四页，共七十三页，编辑于2023年，星期五随机数应具有的基本特性随机数序列应是独立的、互不相关的(uncorrelated)：即序列中的任一子序列应与其它的子序列无关；长的周期(longperiod)：均匀分布的随机数应满足均匀性(Uniformity)：随机数序列应是均匀的、无偏的，即：如果两个子区间的“面积”相等，则落于这两个子区间内的随机数的个数应相等。例如：对[0,1)区间均匀分布的随机数，如果产生了足够多的随机数，而有一半的随机数落于区间[0,0.1]不满足均匀性如果均匀性不满足，则会出现序列中的多组随机数相关的情况均匀性与互不相关的特性是有联系的实际应用中，随机数都是用数学方法计算出来的，这些算法具有周期性，即当序列达到一定长度后会重复；第二十五页，共七十三页，编辑于2023年，星期五有效性（Efficiency):模拟结果可靠模拟产生的样本容量大所需的随机数的数量大随机数的产生必须快速、有效，最好能够进行并行计算。第二十六页，共七十三页，编辑于2023年，星期五随机数与随机数发生器得到一个可能的随机数序列，是在计算机上实现MonteCarlo方法的关键随机数的产生方法：[0,1]区间上均匀分布的随机数是MonteCarlo模拟的基础，服从任意分布的随机数序列可以用[0,1]区间均匀分布的随机数序列作适当的变换或舍选后求得。(0,1)2-1(-1,1)利用随机数表，如Tippett于1972年发表的随机数表；占用太多的计算机内存采用物理方法，如利用电子线路的热噪声等；昂贵而且不便重复:第二十七页，共七十三页，编辑于2023年，星期五伪随机数（Pseudo-RandomNumber)递推到一定次数后，出现周期性的重复现象。利用数学递推公式一旦公式和初值定下来，整个随机数序列便被确定下来，而且每一个随机数只被它前面的那个数唯一确定，因此这类随机数并不是真正的随机数。第二十八页，共七十三页，编辑于2023年，星期五MonteCarlo方法基本思想当所求的问题是某种事件出现的概率，或是某个随机变量的期望值时，它们可以通过某种“随机试验”的方法，得到这种事件出现的频率和概率，或者得到这个随机变量的统计平均值，并用它们作为问题的解。MonteCarlo方法解决的问题问题本身是确定性问题，要求我们去寻找一个随机过程，使该随机过程的统计平均就是所求问题的解。问题本身就是随机过程，我们可以根据问题本身的实际物理过程来进行计算机模拟和跟踪，并采用统计方法求得问题的解。第二十九页，共七十三页，编辑于2023年，星期五MonteCarlo方法的特点计算的收敛性和收敛速度均与问题的维数无关，适合解决高维问题。对问题的适应能力强。收敛速度仅为样本数的-1/2次，因而计算耗时大。第三十页，共七十三页，编辑于2023年，星期五MonteCarlo方法的应用举例：计算积分：常用的积分方法求解：将积分区域[a,b]均匀地划分成N各分区间，则积分结果可近似地表示成：Δx=(b-a)/N第三十一页，共七十三页，编辑于2023年，星期五简单的MonteCarlo积分方法求解：利用均匀分布的随机数发生器，从[a,b]区间产生一系列随机数xi，i=1,

2,

...,

N其中X为均匀分布，并且X[a,b]近似求解E[g(X)]:近似求解积分:随机抽样第三十二页，共七十三页，编辑于2023年，星期五当我们用简单MonteCarlo计算积分时，若该函数为常数函数，g(x)=constant，则取样数不管多少，准确度为100％。如果在积分区间内，g(x)为一平滑函数，则简单MonteCarlo方法较为准确，反之，如果g(x)的变动很剧烈，则简单MonteCarlo方法的误差会变大。说明：第三十三页，共七十三页，编辑于2023年，星期五重要性MonteCarlo抽样方法在g(x)变化剧烈时，如果以MonteCarlo方法取样，最好依据g(x)的大小来决定取样率。当|g(x)|的值较大时，对∫g(x)dx的贡献也较大，如果没被选中，则结果的误差极大。解决方式：改变

x被选中的机率，让|g(x)|值较大的点被选中的机率增加。采用权重分布函数(Weightdistributionfunction)w(x):决定每个x被选中的机率。重要性抽样的定义：根据一定的分布形式进行的随机抽样。第三十四页，共七十三页，编辑于2023年，星期五w(x)必须归一化，即在积分区间内∫w(x)dx=1。由于

x的选取已被w(x)扭曲，所以计算积分时要把这部分[还]回去：若一共取样了N个x，则积分值为：重要性MonteCarlo抽样方法第三十五页，共七十三页，编辑于2023年，星期五

MetropolisMonteCarlo方法我们所模拟的系统最终要达到的平衡分布是Boltzman分布：Boltzmann概率分布函数：我们如果能够产生这种分布，我们就能够计算系统的大多数性质，但这是不可能的，因为我们不知道Z的值，但是对于任意两个状态，我们有：可以在相空间中构造一个马尔科夫链，使相空间中的样本点随着链的增长逐步趋近于Boltzman分布。第三十六页，共七十三页，编辑于2023年，星期五一个序列x0,x1,x2,…,xn,如果对任何n都有：则此序列是一个Markov链。要求：任何一次实验的结果依赖于前一次的试验，并且近依赖于前一次的试验。马尔科夫(Markov)链可以证明：通过构造Markov链，体系中最终的平衡分布就是Boltzman分布第三十七页，共七十三页，编辑于2023年，星期五Metropolis平衡条件(Detailedbalancecondition)：平衡条件：系统处于状态X的概率正比于其Boltzman因子：如果是对称的:第三十八页，共七十三页，编辑于2023年，星期五

MetropolisMonteCarlo方法的算法：给出一个初始状态，并计算系统的能量Eold随机产生一个新状态，并计算新系统的能量Enew如果ΔE＝(EnewEold)<0,则接受新状态并回到b)如果ΔE＝(EnewEold)>0,则计算Boltzman因子：在(0,1)区间上产生一个均匀分布的随机数；如果则接受新状态并回到b)否则保留原值并回到b)第三十九页，共七十三页，编辑于2023年，星期五1.正则系综蒙特卡罗模拟方法(CanonicalMCSimulation)具有确定的粒子数N、温度T和体积V第四十页，共七十三页，编辑于2023年，星期五对于含N个粒子的系统，位型（构型）的配分函数：某个特定构型的发生概率为PNVT(rN)1.正则系综蒙特卡罗模拟方法(CanonicalMCSimulation)第四十一页，共七十三页，编辑于2023年，星期五

MonteCarlo模拟中任一物理量的计算：位型积分概率密度系统处于位型{rN}的概率密度第四十二页，共七十三页，编辑于2023年，星期五给出一个初始状态，并计算系统的能量Uold随机产生一个新状态，并计算新系统的能量Unew如果ΔU＝(UnewUold)<0,则接受新状态并回到b)如果ΔU＝(UnewUold)>0,则计算Boltzman因子：在(0,1)区间上产生一个均匀分布的随机数

；如果则接受新状态并回到b)否则保留原值并回到b)正则系综MC模拟算法的组织：第四十三页，共七十三页，编辑于2023年，星期五正则系综MC模拟算法的流程图：给定每个分子的初始位置，ri(0)随机选取一个分子，并随机移动到新的位置计算移动前后的系统能量变化ΔU拒绝移动ΔU0?Exp(ΔU)(0,1)?统计系统的热力学性质及其它物理量统计性质不变？

打印结果，结束YesNoNoNo接受移动YesYes大约循环107到108次第四十四页，共七十三页，编辑于2023年，星期五MonteCarlo模拟中几个热力学量的计算：N个粒子系统中的总势能：假设采用截断势能函数：Uc：截断范围内的总势能；Ulrc：截断半径外对势能的长程校正(Long-rangecorrection)对于LJ流体：第四十五页，共七十三页，编辑于2023年，星期五

含N个粒子系统中的压力：Wc：截断范围内的总维里项(Virial)；Plrc：截断半径外对压力的长程校正第四十六页，共七十三页，编辑于2023年，星期五ReadsimulationparametersStartInitializepositionsofallparticlesNewsimulation?ReadoldconfigurationMonteCarloloopStopyesnoMonteCarloloopSubroutineStartStopTrialmoveSatisfy

Metropolisrule?AcceptthetrialmoveUpdateenergyandvirialSamplethepressureEndofsimulation?yesnoyesnoMainprogram正则系综MC模拟程序基本结构：第四十七页，共七十三页，编辑于2023年，星期五正则系综MC模拟程序F11讲解（LJ,NVT）：**READINPUTDATA**初始状态：READ(*,‘(A)’)Title!运行作业题目READ(*,*)NStep!运行步数READ(*,*)Iprint!打印步数READ(*,*)Isave!保存步数READ(*,*)Iratio!调整步数READ(*,‘(A)’)CNFile!位型文件READ(*,*)Dens!对比密度READ(*,*)Temp!对比温度

READ(*,*)Rcut!对比截断半径第四十八页，共七十三页，编辑于2023年，星期五无因次量：第四十九页，共七十三页，编辑于2023年，星期五正则系综MC模拟程序F11讲解（LJ,NVT）：量纲变换：Beta=1.0/TempSigma=(Dens/Real(N))**(1.0/3.0)Rmin=0.70*Sigma！判断粒子发生重叠时的距离Rcut=Rcut*Sigma！截断半径DRmax=0.15*Sigma！随机移动的最大距离DensLJ=DensDens=Dens/(Sigma**3)IF(Rcut.GT.0.5)Stop'Cut-OffTooLarge'模拟盒子的边长为1第五十页，共七十三页，编辑于2023年，星期五**ReadInitialConfiguration*CallReadCN（CNFile）正则系综MC模拟程序F11讲解（LJ,NVT）：初始位型：参阅程序F23CallFCC需要自己给定所有粒子初始位置面心立方(face-centeredcubic,FCC)：第五十一页，共七十三页，编辑于2023年，星期五正则系综MC模拟程序讲解（LJ,NVT）：长程校正：Sr3=(Sigma/Rcut)**3Sr9=Sr3**3Vlrc12=8.0*Pi*DensLJ*Real(N)*Sr9/9.0Vlrc6=-8.0*Pi*DensLJ*Real(N)*Sr3/3.0Vlrc=Vlrc12+Vlrc6Wlrc12=4.0*Vlrc12Wlrc6=2.0*Vlrc6Wlrc=Wlrc12+Wlrc6第五十二页，共七十三页，编辑于2023年，星期五算法：能量求和：CallSumup(Rcut,Rmin,Sigma,Ovrlap,V,W)If(Ovrlap)Stop'OverlapInInitialConfiguration'Vs=(V+Vlrc)/Real(N)Ws=(W+Wlrc)/Real(N)Ps=Dens*Temp+W+WlrcPs=Ps*Sigma**3第五十三页，共七十三页，编辑于2023年，星期五**CheckForAcceptance**DeltV=Vnew-VoldDeltW=Wnew-WoldDeltVb=Beta*DeltVIf(DeltVb.Lt.75.0)ThenIf(DeltV.Le.0.0)ThenV=V+DeltVW=W+DeltWRX(i)=RXInew…………..Acatma=Acatma+1.0Elseif(Exp(-DeltVb).Gt.Ranf(Dummy))ThenV=V+DeltVW=W+DeltWRX(i)=RXInew……………Acatma=Acatma+1.0EndifEndifAcm=Acm+1.0算法：Metropolis算法第五十四页，共七十三页，编辑于2023年，星期五算法：势能的计算DO100J=1,NIF(I.NE.J)ThenRXIJ=RXI-RX(J)………RXIJ=RXIJ-Anint(RXIJ)……………RIJSQ=RXIJ*RXIJ+RYIJ*RYIJ+RZIJ*RZIJIF(RIJSQ.LT.RcutSQ)THENSR2=SIGSQ/RIJSQSR6=SR2*SR2*SR2VIJ=SR6*(SR6-1.0)WIJ=SR6*(SR6-0.5)V=V+VIJW=W+WIJENDIFENDIF100ContinueV=4.0*VW=48.0*W/3.0最小影像原理第五十五页，共七十三页，编辑于2023年，星期五算法：Metropolis算法RXIOld=RX(I)…………..**CalculateTheEnergyOfIInTheOldConfiguration**CallEnergy(RXIold,RYIold,RZIold,I,Rcut,Sigma,Vold,Wold)C**MoveIAndPickupTheCentralImage**RXInew=RXIold+(2.0*Ranf(DUMMY)-1.0)*DRmax………………..RXInew=RXInew-Anint(RXInew)………………**CalculateTheEnergyOfIInTheNewConfiguration**CallEnergy(RXInew,RYInew,RZInew,I,Rcut,Sigma,Vnew,Wnew)随机移动第五十六页，共七十三页，编辑于2023年，星期五**AdjustMaximumDisplacement**Ratio=Acatma/Real(N*Iratio)If(Ratio.Gt.0.5)ThenDRmax=DRmax*1.05ElseDRmax=DRmax*0.95EndifAcatma=0.0EndifIf(Mod(Step,Iprint).Eq.0)Then算法：步长的调整第五十七页，共七十三页，编辑于2023年，星期五分子模拟中径向分布函数g(r)的计算：计算表达式：物理意义：流体中距一个分子为r处出现另一个分子的几率密度，它反映流体中短程有序的特点。rr+r在模拟中通常在盒子长度一半的范围内，考察g(r)随距离的变化。第五十八页，共七十三页，编辑于2023年，星期五分子模拟中径向分布函数g(r)的算法：第一步：计算球壳层间的距离，并初始化一些变量：Ngr:g(r)的统计次数Delg:球壳层间的距离nhist:球壳层的层数模拟开始时需要给定：nhist以及统计g(r)的步数间隔。第五十九页，共七十三页，编辑于2023年，星期五分子模拟中径向分布函数g(r)的算法：第二步：统计g(r):第六十页，共七十三页，编辑于2023年，星期五分子模拟中径向分布函数g(r)的算法：第三步：系综统计平均计算g(r):第六十一页，共七十三页，编辑于2023年，星期五第六十二页，共七十三页，编辑于2023年，星期五2.巨正则系综蒙特卡罗模拟方法(GrandCanonicalMCSimulation)恒定V,T,和

，体系的粒子数发生波动；可用于预测