中国矿业大学分子模拟课件第十章

分子模拟牛继南[e-mail]2011.3第十章分子动力计算结果分析建模获得力场参数计算分析1/21/210.1运动轨迹分析•三种方式：1：直接动画显示正戊烷2：随时间的变化曲线θ和ϕ的协同运动3：统计分布曲线演示相邻两个C-C-C角的分布状况θ1θ2θ1和θ2的统计分布1001051101151201251300.000.020.040.060.080.10ProbabilityAngle(degree)θ2θ110.2热力学特性计算1：系统的温度可以根据系统的平均动能得到系统的平均温度。其中：K为动能T为温度f为系统的自由度34036038040042044046048050001234567891011121314151617181920Temperature(K)Time(ps)ForciteAnalysis-TemperatureProfileRunningaverage2：动能的统计涨落系统动能的统计涨落(fluctuation)定义为：通过统计涨落可以计算出系统的定容比热容（heatcapacityatconstantvolume）Cv：式中N为原子总数，K为系统动能，kB为玻尔兹曼常数。-100010020030001234567891011121314151617181920Energy(kcal/mol)Time(ps)ForciteDynamics-EnergiesPotentialenergyKineticenergyNon-bondenergyTotalenergy3：压力和应力张量压力为单位面积上所受到的力，如标准气压为1.013kPa。实际中系统各个方向上所受的力并不相同，通常用一个张量表示：对于均匀液体或气体，压力在各个方向相同，张量为：P为静压。应力张量•常用的压力单位为GPa，定义正的压力值为压缩系统向内的力所致，定义正的应力值为拉伸系统向外的力所致，两者采用的标准相反。•系统的压力可以由体积、温度和位力（virial）计算：•其中位力W为：或如果将瞬时温度带入：得：因为压力为一张量，故可写为：•可以按照这两个张量求出压力张量，通常的所说静压力为压力张量迹的和的三分之一。•对于某一种材料，可以通过不同压力下的分子动力学计算获得体积或晶胞参数在不同压力或应力的改变量，从而建立应力-应变曲线（strain-stresscurve）。•因为分子动力学计算当中通常会采用周期性边界条件，计算系统的位力时，必须考虑原系统和镜像系统的作用，所以i≠j，即：-1.2-1.0-0.8-0.6-0.4-0.20.00.20.40.60.81.001234567891011121314151617181920Pressure(GPa)Time(ps)ForciteAnalysis-PressureProfileRunningaverage10.3径向分布函数计算•系统中参考分子周围r~r+dr范围内的分子数目为dN，则定义径向分布函数g(r)为：•其中ρ为系统的密度，若系统的总分子数为N，可得：•因为：•所以径向分布函数可以解释为系统的区域密度和平均密度之比。•参考分子附近的区域密度不同于系统密度，但与参考分子距离远时区域密度应与平均密度相同，即当r值大时径向分布函数的值应当接近于1.•分子动力计算径向分布函数的方法为：024681012141618202224262801234567891011121314g(r)r(Angstrom)ForciteAnalysis-RDF(C)Total01234501234567891011121314g(r)r(Angstrom)ForciteAnalysis-RDF(H)Total•径向分布函数反映了液体或气体分子的聚集特性，可借此了解流体的结构。将g(r)在0~rmax范围内积分，可以得到此范围内的分子平均数目。如水分子周围的水和数目即可利用这种方法求出。10.4与时间有关的物理量计算•分子动力计算除了计算系统的平均值以外，最重要的是计算系统的各种动态特性。1：均方位移分子动力计算当中，分子由起始位置不停的移动，每一瞬间的位置都不相同，如果以表示时间t时粒子i的位置，则粒子位移平方的平均值称为均方位移（meansquaredisplacement，MSD），即：•根据统计原理，只要分子数目足够多，计算的时间足够长，系统的任何一个瞬间都可以作为起点，计算的平均值相等，因此计算MSD时常采用多时间起点：设一共运行n步，各步位移，将此轨迹分为相等数目的两部分:计算MSD时每次均取n/2组数据的平均:01234567891011121314150.00.10.20.30.40.50.60.70.80.91.01.11.21.31.41.51.61.71.81.92.0MSD(Angstrom^2)Time(ps)ForciteAnalysis-MeanSquareDisplacement(C)TotalMSD024681012141618200.00.10.20.30.40.50.60.70.80.91.01.11.21.31.41.51.61.71.81.92.0MSD(Angstrom^2)Time(ps)ForciteAnalysis-MeanSquareDisplacement(H)TotalMSD•根据爱因斯坦扩散定律，有：•式中D为粒子的扩散系数，因此当模拟的时间很长时，均方位移对时间曲线的斜率即为6D，由此可求出粒子的扩散速率。DC=1.03Å2/sDH=1.15Å2/s2：相关函数设A，B为系统两个不同的物理量（位置、速度等），则A或B的自身相关函数(autocorrelationfunction)可表示为：相关函数的物理意义为一个物理量随时间改变后与其起始时的相关性。因为物理量的平均值不随选择时间起点而改变，所以上式中T可为任意时刻。•A和B的交互相关函数（crosscorrelationfunction）的表达式为：•相关函数是表示物理量与物理量间与时间的相关性。通常定义归一化的自身相关函数为：•由分子动力模拟所储存的轨迹计算相关函数的方法与计算均方位移一样，应将轨迹的格点均视为时间的起点。•设分子动力计算系统含有N个分子，共收集了n步轨迹，则：•下图所示为存储轨迹中t=4的相关性。•设积分的步长为δt，则经过τδt时的相关函数为：•τmax为选取的最大轨迹桢差。•对于相关函数的计算，通常有两种方法：1.直接计算法缺点：占用的存储空间大；2.快速傅里叶变换法复杂，但更有效FFT(FastFourierTransform)•如果不考虑归一化系数，则时间相关函数可改写为：•'表示未含归一化常数，2τrun表示计算了2倍所储存的资料点。实际操作中，可在一组资料点后加上一组零，以满足此要求。上式还可写为：•根据FFT的卷积(convolution)定理，若：•则：•式子中的表示时间空间向频率空间的转换τ→ν，即：•如果从ν→τ，则：•根据这个关系，相关函数可以转化为：•所以：•因此FFT计算时间相关函数的步骤为：1.将τrun个MD资料加上τrun个零；2.利用FFT方法求出A(τ)的变换A(ν)；3.计算；4.利用反傅里叶变换计算5.乘上归一化常数，求得。-1001020304050607080900.00.10.20.30.40.50.60.70.80.91.01.11.21.31.41.51.61.71.81.92.0VACF(Angstrom^2/ps^2)Time(ps)ForciteAnalysis-VelocityAutocorrelationFunction(C)VACF短时间内，粒子速度和其初速度关联性强；当τ足够大时，粒子的速度已经和其初速度没有关联性-300-200-100010020030040050060070080090010000.00.10.20.30.40.50.60.70.80.91.01.11.21.31.41.51.61.71.81.92.0VACF(Angstrom^2/ps^2)Time(ps)ForciteAnalysis-VelocityAutocorrelationFunction(H)VACF短时间内，粒子速度和其初速度关联性强；当τ足够大时，粒子的速度已经和其初速度没有关联性•速度自相关函数可以计算粒子的扩散系数，关系为：•对C和H的VACF积分，得到：•DC=0.969Å2/s•DH=1.139Å2/s•表明在系统中H原子运动的速度要快于C原子。DC=1.03Å2/sDH=1.15Å2/s由MSD求出的为速度相关函数除了计算粒子的扩散速率以外，还有一些重要的应用。•与粒子运动相关的运动模式有其对应的频率，它们通常在红外光谱范围（400~4000cm-1）。与原子运动相关的频率可以由速度相关函数的傅里叶变换得到：•其中ω为频率，Re表示实数部分，I(ω)为频率的能谱密度,根据轨迹综合上式可写为：010203040010020030040050060070080090010001100120013001400150016001700Intensity(Angstrom^4/ps^2)Frequency(1/cm)ForciteAnalysis-VACFPowerSpectrum(C)Total0100200300400500600700800010020030040050060070080090010001100120013001400150016001700Intensity(Angstrom^4/ps^2)Frequency(1/cm)ForciteAnalysis-VACFPowerSpectrum(H)Total1.阻碍运动2.可以通过振动峰位置，理解分子间作用。•分子偶极的相关函数可以用来计算红外光谱，光谱的线状函数为：•其中为分子的偶极矩。•红外光谱仪检测到的频率位置和光谱计算得到光谱位置应该相同，因此可以对比验证所用力场参数的正确性。0.000000.000020.000040.000060.000080.000100.000120.000140.000160.000180.00020010020030040050060070080090010001100120013001400150016001700Intensity(ps^2)Frequency(1/cm)ForciteAnalysis-DACFPowerspectrum(all)Total•可以计算拉曼散射光谱，拉曼光谱的函数为：•其中为归一化的转动相关函数（rotationalcorrelationfunction）-0.10.00.10.20.30.40.50.60.70.80.91.00.00.10.20.30.40.50.60.70.80.91.01.11.21.31.41.51.61.71.81.92.0RTCFTime(ps)ForciteAnalysis-RotationalTimeCorrelationFunction(dis)M1M2质心到C3原子的距离m1和m2是两种不同的转动相关函数•还可计算粘度和热导率等•粘度：其中•热导率：•其中这两个值的误差通常较大，是因为模拟的体系太小的缘故