清华《分子动力学模拟入门》

多体系统—分子动力学Multi-bodySystemMolecularDynamics引言物质基本构成—分子、原子•在分子、原子这个微观水平上来考察物质：多体世界•查清楚微观世界，宏观就清楚了从微观考虑问题的现实可行性从微观考虑问题的必要性•物性的观测性参数：热传导、温度、压力、粘性、...…微观处理的前提•已知微观粒子间的相互作用假设•分子为球，惰性，分子间的作用只取决于分子间的距离分子动力学(MolecularDynamics，MD)MD的应用领域：物理、化学、生物、材料等MD方法能实时将分子的动态行为显示到计算机屏幕上,便于直观了解体系在一定条件下的演变过程MD含温度与时间,因此还可得到如材料的玻璃化转变温度、热容、晶体结晶过程、输送过程、膨胀过程、动态弛豫(relax)以及体系在外场作用下的变化过程等水和离子在微小硅孔中的运动聚乙烯的结晶MD的基本原理用牛顿经典力学计算许多分子在相空间中的轨迹•求解系统中的分子或原子间作用势能和系统外加约束共同作用的分子或原子的牛顿方程。•模拟系统随时间推进的微观过程。•通过统计方法得到系统的平衡参数或输运性质•计算程序较为复杂，占用较多内存MD的主要步骤选取要研究的系统及其边界，选取系统内粒子间的作用势能模型设定系统中粒子的初始位置和初始动量建立模拟算法，计算粒子间作用力及各粒子的速度和位置当体系达到平衡后，依据相关的统计公式，获得各宏观参数和输运性质分子间势能及相互作用N个粒子系统的总势能NjiijNNNNrVrVrVrVrVrVrVV11,-22311312)()()()()()()(dVrVr刚球模型()0rdVrrd()()dVrr斥力力心点模型drV-()()rdVrdrdrSoutherland模型分子间势能及相互作用Lennard-Jones势能能量尺度；长度尺度为方便，时常归一化：记；6124)(rrrV612/1/14)(rrrVVV/rr/612114)(rrrVV(r)F(r)排斥力V(r),F(r)r吸引力分子间势能及相互作用一些气体的参数kB=1.38x1023(J/K):Boltzmann常数NeonArgonKryponXenonNitrogen(nm)0.2750.34050.3600.4100.370/kB(K)36119.817122195分子间势能及相互作用相互作用•标量形式：•直角坐标：至此，各粒子间相互作用已知，可进行模拟了)(rVrmf()(),yxffVrxVrymrrmrr)(rVmF模拟的数学方法Euler法和Euler-Cromer方法？•不能用：不能保持总能量守恒•Verlet算法：速度形式12121212111212112nnnnnnnnnnnnnnxxvtatvvaatvvatvvat模拟的数学方法Leap-frog算法：x的截断误差为，v的截断误差1/21/211/21/22nnnnnntnnnvvtaxxtvvva4()Ot2()Ot模拟的数学方法边界条件•模拟能力限制，不能模拟大量分子，只能模拟有限空间中的有限个分子：有限空间边界•固体（刚性）边界条件不仅仅有分子间的相互作用，还引入了壁面的作用分子量大时，壁面作用可忽略不计233611aaaN和壁面作用分子数壁面积总分子数体积3aN总分子数壁面积和壁面作用分子数模拟的数学方法取~，前比值为0.2~0.01。取前值，模拟粗糙；取后值，模拟计算量太大处理方法：使用周期性边界条件•周期性边界条件210N610a模拟的数学方法两个不同粒子在x或y方向上的最大分离距离为a/2•最小像约定：两粒子分离距离最大分离距离，相互作用力可以忽略，而加入其中像粒子之一相互作用力来考虑?2a2/aa模拟的数学方法考虑的粒子总数不变初始条件•随机初始条件给法之一要求大小：条件一：规则给法条件二：随机给法max||Vv)1random2(maxVv模拟的数学方法random：随机数产生函数，产生(0，1)之间的随机数。方向(按球坐标给法)：分量：xyz))5.0random(2arccos()]5.0random(2[arccos)0.5random(sign180180;1800cossinvvxsinsinvvycosvvz模拟微观量•温度根据统计热力学，平衡态下经典系统的能量中的每一个二次项具有平均值kBT/2,即注意：上式在系统质心速度为0时适用iiiBvmTNkd2212空间维数粒子个数：取时间平均模拟问题：如何给定系统的初始条件，得到所需要的平衡态温度Teq？解决方法之一：速度标定法任给初始条件，模拟到平衡，得到系统平衡态温度T。一般TTeq。令用速度再模拟直到平衡，若所得温度仍不等于Teq，再进行上述过程TTfeq/1/diifvv给定初始条件：xi，vi计算到平衡态|T-Teq|≤εf=Teq/Tvi=vif1/d计算结束YN其它方法：Gaussian热浴法（约束温度调节方法）其基本原理在运动方程中加入“摩擦力”项，并将其与粒子速度联系起来。平衡态时，系统温度不变，因此dEk/dt=0iiifmamv0iiiva2iiiiifvmv宏观性质的统计系统的势能系统的内能系统的总能E=Ep+Ek系统的温度1()pijijNEVr22kiipEm21iiiBTmvdNk模拟•热容定义热容计算系统在温度T和T+T时的总能ET、ET+T，VvTECE：系统总能TEECTTTv模拟•压强对壁面的压强t时间里作用在单位面积壁上的压力tiv时刻，速度为ttiv时刻，速度为iiiitvvmdAp1刚性壁模拟•粒子速度分布选速度间隔v，模拟nt个时间步，记录在每个速度间隔中的粒子数，最后归一化。vv)(vNvvvv2模拟气、液状态方程维里定理(VirialTheorem)iiiBFrdTNkpV1压强体积粒子i的位置矢量粒子i所受到的其它粒子的合相互作用力温度的模拟可得此项在温度的模拟基础上再模拟此项模拟例：用此可确定高密度气体和液体状态方程(vanderWaals方程)中的系数理想气体状态方程在高密度情况下不可用确定系数a和b21abTkpB气体密度CASE(1)-CouetteFlowSizeofdomainis:12.51x7.22x16.711,1,1.5,0.6,0.02,1.4ll6126124()ijijijijErr6127124(612)ijijijijdEFdrrrIfδislessthan0,thenthetwospeciesareimmiscible.i,jrepresentdifferentspeciesCASE(1)-CouetteFlowztimeCASE(1)-CouetteFlowzCASE(2)-ContactAngleSimulationMass:m[1]=1,m[2]=8,m[3]=0.8L=25.05,W=6.56H=10.29T=1.2110.21110.211ij126()ijijLJijijVrrri,j=1,2,3,1–redfluid,2--wall,3–greenfluidVCASE(3)-Rayleigh-Taylorinstabilitygeneration(重力场)Dzwinel,W.,Alda,W.,Pogoda,M.,andYuen,D.A.,2000,Turbulentmixinginthemicroscale:a2Dmoleculardynamicssimulation,PhysicaD,Vol.137,pp.157-171.Fig.1.ThesnapshotsofMDsimulationsoftheR–Tinstabilityfortwoparticlesystems:(A)closed;(B)open(theheavyfluidiscolouredinlightgreywhilethelighteroneisdarkgrey,theblackpartofthefigureisempty,i.e.,itdoesnotcontainparticles).Fig.2.TheeffectofconverginggeometryobtainedbyMDsimulationofonemillionparticlesinthemicroscale.0.5μmCASE(4)Typicaltranslocationevent1.A1.4Vbiasappliedtomembrane.2.20base-pairfragmentofdoublestrandedDNAplacedinfrontofananopore.3.EndofDNAnearesttotheporeispulledintotheporebyitschargedbackbone(a,b)4.Systemreachesameta-stablestate(c)andtranslocationhalts.5.Base-pairsstarttosplit.Somefreednucleotidesadheretoporesurface.6.Voltageincreasedmomentarilytodrivesystemoutofmetastablestate.7.DNAexitspore.Oneofthebasesholdsfirmlytotheporesurface.8.After50ns,mostofDNAhasleftpore.Nineoftwentybasepairsaresplit.BubblenucleationonsolidsurfacesMaruyama,S,andKimura,T.,2000,AMolecularDynamicsSimulationofaBubbleNucleationonSolidSurface,HeatandTechnology,Vol.18,pp.69-73.MoleculartransportindropletsMaruyama,S.,Matsumoto,S.,andOgita,A.,1994,SurfacePhenomenaofMolecularClustersbyMolecularDynamicsMethod,ThermalScienceandEngineering,Vol.2,No.1.FullViewSlicedViewVoidView参考书目D.C.Rapaport.TheArtofMolecularDyna-micsSimulation.LAMMPS开源的MD模拟软件版本•Fortran77＋MPI•Fortran90＋MPI•C++＋MPI可单机运行，也可并行计算可以模拟气体、液体、固体SamplesMDProgramcrystal.f90•生成粒子•生成数据文件md1.f90•MD模拟