分子动力学课程小结.docx

《分子动力学课程小结.docx》由会员分享，可在线阅读，更多相关《分子动力学课程小结.docx（24页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、分子动力学课程小结分子动力学模拟课程小结一分子动力学的基本原理在分子动力学模拟中，体系原子的一系列位移是通过对牛顿运动方程积分得到的，结果是一条运动轨迹，它表明了系统内原子的位置与速度怎样随时间而发生变化。通过解牛顿第二定律的微分方程，能够获得原子的运动轨迹。方程如下：这个方程描绘了质量为mi的原子i在力Fi的作用下，位置矢量为ri时的运动方程。其中，Fi能够由势函数U的梯度给出：系统的温度则与系统中全部原子的总动能K通过下式相联络：N是原子数，Nc是限制条件，kB是波尔兹曼常数。二.MD模拟的积分算法为了得到原子的运动轨迹，能够采用有限差分法来求解运动方程。有限差分法的基本思想就是将积分分成

2、很多小步，每一小步的时间固定为t。用有限差分解运动方程有很多方法，所有的算法都假定位置与动态性质速度、加速度等能够用Taylor级数展开来近似：在分子动力学模拟中，常用的有下面的几中算法：1.Verlet算法运用t时刻的位置和速度及tt时刻的位置，计算出tt时刻的位置：两式相加并忽略高阶项，能够得到：速度能够通过下面方法得到：用tt时刻与tt时刻的位置差除以2t：同理，半时间步t+t时刻的速度可以以算：Verlet算法执行简单明了，存储要求适度，但缺点是位置r(tt)要通过小项与非常大的两项2r(t)与r(tt)的差相加得到，容易造成精度损失。另外，其方程式中没有显示速度项，在没有得到下一步的

3、位置前速度项难以得到。它不是一个自启动算法：新位置必须由t时刻与前一时刻tt的位置得到。在t0时刻，只要一组位置，所以必须通过其它方法得到tt的位置。一般用Taylor级数：2.VelocityVerlet算法3.Leap-frog算法为了执行Leap-frog算法，必须首先由t-0.5t时刻的速度与t时刻的加速度计算出速度v(t+t)，然后由方程计算出位置r(t+t)。T时刻的速度能够由：得到。速度蛙跳过此t时刻的位置而得到t0.5t时刻的速度值，而位置跳过速度值给出了tt时刻的位置值，为计算t0.5t时刻的速的作准备，依此类推。其缺点是位置与速度不同步。这意味着在位置一定时，不可能同时计算

4、动能对总能量的奉献。三.分子动力学计算的时间间隔时间间隔t在积分算法中是一个非常重要的参数。为了充分利用CPU时间，尽量选择比拟大的时间间隔，但是假如时间间隔太大，就会造成积分经过的不稳定性和不准确性。时间间隔的设置同时依靠于算法和模型的情况。模型本身给时间间隔带来的最大的限制就是最高频率的运动。由于Verlet算法要求在每个时间间隔内模型的速度和加速度保持一边，时间间隔就应该低于振动周期的八分之一到特别之一。对大多数的有机模型来讲，最高的振动频率是C-H键的伸缩振动，其振动周期的数量级为1014s。这样，时间间隔就应该是0.5-1fs左右。假如采用受约束的SHAKE或者RATTLE算法，能够

5、使用更长的时间间隔。假如研究对象是液态或者固态简单模型，对体系内作用形式不感兴趣，可以以采用一些更长的时间间隔，比方20fs。对离子态的材料模型，5fs左右是适宜的。时间间隔必须跟选择的算法相匹配。比方，ABM4算法的时间间隔应该是Verlet算法的一半左Runge-Kutta-4算法则需要比其他算法更短的时间间隔。四.生物大分子的互相作用势函数生物大分子势一个自有度很大，构造复杂的体系，各种物理经过极为丰富。其中的互相作用规律还难以进行统一的处理。但其物理图像是明晰的，电子的运动和骨架有效的分开，成键电子是定域的，价键的概念仍然有效。荷电集团之间存在静电互相作用。原子或集团正负电荷中心不重合

6、构成电多极矩并介入静电互相作用。原子或集团的电子云受静电互相的诱导构成电诱导偶极由此产生色散力。静电、范德华互相作用称非键互相作用。此外有共价键德振动及单键德旋转绕共价双键德扭曲等称为成键互相作用。典型德生物分子动力学的模型势函数能够写成一下形式：第一项求和是共价键的键能，第二项求和是键角能，第三项求和是非正常二面角能量，第四项求和是正常二面角能量。Kb，k,k，k代表式中前四项成键互相作用的力常数。式中最后一项式体系中的非键互相作用项，(i,j)表示对非键对求和，C12、C6式范德华互相作用势中的系数。势函数中参数qi，qj分别为第i荷第j个原子的电荷。但是，对于疏水互相作用，溶剂的极化荷对

7、蛋白质内部静电屏蔽，都需要对势函数加以修正和补充。若要了解溶液中分子的动力学性质则须将模拟的体系放到充分多的水分子做充分动力学模拟。五.分子动力学模拟的启动进行MD模拟，必须首先建立系统的初始构造。初始构造能够通过实验数据、理论模型、或两者的结合来获得。除此之外，还能够给每个原子赋予初速度，它能够从一定温度下的Maxwell-Boltzmann分布来任意选取：Maxwell-Boltzmann分布给出了质量为mi的原子i在温度T下沿x方向速度为vix的概率。Maxwell-Boltzmann分布势一中Guassian分布，它能够用随机数发生器得到。大多数随机发生器产生的随机数均匀分布在01之间

8、，但是能够通过变换得到Guassian分布。均值为和波动值为2的Guassian分布的概率为：一种方法势首先产生两个在01之间的随机数1和2。运用下列式子可产生两个数下x1、x2：另一种方法势先产生12个随机数1、2、12，然后计算：这两种方法产生的随机值都服从均值为零，偏差为一个单位的正态分布。初始速度经常被校正以知足总动量为零，为了是总动量为零，分别计算沿三方向的动量总和，然后用每一方向的总定量除以总质量，得到一速度值。用每个原子的速的减去此速度值，即可保持系统的总动量为零。在建立了系统的初始位形和赋予初始速度后，就具备分子动力学模拟的初步条件了。在每一步中，原子所收到的力通过对势函数的微

9、分能够得到。然后根据牛顿第二定律，计算加速度，再由以上提供的算法即可进行连续的模拟计算了。六.分子动力学模拟的系综采用MD模拟，必须再一定的系综下进行，经常用到的系综包括微正则系综、正则系综、等温等压系综和等温等焓系综。1.微正则系综(NVE)是孤立的、保守的系统的系综，在这种系综中，系统沿着相空间中的恒定能量轨道演化。在分子动力学模拟的经过中，系统中的原子数(N).体积V、和能量E都保持不变。2.正则系综NVTNVT保持不变，并且总动量为零。恒温下，系统的总能量不是一个守恒量，系统要与外界发生能量交换。保持系统的温度不变，通常运用的方法是让系统与外界的热浴处于热平衡状态。由于温度与系统的动能

10、有直接的关系，通常的做法是把系统的动能固定在一个给定值上。3.等温等压系综NPTNPT保持不变，这种系综是我们常见的系综，很多分子动力学模拟都要在从系综下进行。这时，要保证系统的温度恒定，还要保持它的压力恒定。温度恒定和以前一样，是通过调节系统的速度或加一约束力来实现的。而对压力进行调节，就比拟复杂。由于系统的压力P与其体积V是共轭量，要调节压力值能够通过改变系统的体积来实现。4.等压等焓系综NPHNPH保持不变，焓值通过H=E+P得到。在系综下进行模拟时要保持压力和焓值为固定值。七.边界条件正确处理边界和边界效应对模拟方法时至关重要的，由于它时从模拟相对较少的原子来计算物质的宏观性质的。为了

11、减小有限尺寸的影响，经典的方法是使用周期性边界条件。周期性边界条件使得能够用相对少的粒子数目来真实的模拟大块的体系，使得粒子仿fo处在一个完好的体系中，假设一个用来模拟的立方单胞，使这个单胞在各个方向上都不断重复，看起来象有周期一样。在二维图像中，每一个单胞都被其他的8个单胞所包围；在三维方向上，每一个单胞就会被26个单胞所包围。所有的单胞中的粒子的坐标都能够通过一个整数而得到。当模拟的单胞中的一个粒子由于力的作用而离开这个单胞的时候，九会有另一个和它对应的粒子运动到这个单胞中来，这样，模拟的整个体系的粒子数就会保持不变。在分子动力学模拟中经常采用的边界条件有：1.矩形盒子周期性边界条件2.单

12、斜盒子周期性边界条件3.去头八面体合租周期性边界条件八.MD模拟结果分析方法MD模拟的时间能够到达几百皮秒或纳秒，甚至更长，在运行MD模拟时，体系的速度和坐标被保存下来。在分析时能够对热力学参数进行计算。热力学参数随时间演化的特性能够用图形显示，一个坐标对应时间，另外一个坐标时感兴趣的物理量，如能量、均方根差、原子位置的涨落，此外还能计算出平均构造与实验数据相对照，这些都有助于在原子水平上直观地理解构象地变化。1.均方根差(RMSD)一种广泛用于检验MD模拟正确性地方法时求相对于蛋白质晶体结果地均方根差。目前已知一部分蛋白质地晶体构造与在溶液中的结果存在显著的差异，但人们普通以为这种差异对于大

13、多数蛋白质时非常小的。蛋白质模拟中常有的质量控制时获得既小又稳定的蛋白质骨架原子的RMSD值(通常其中，e(r,R)是把核坐标作为参变量的电子波函数冲N(R)是原子核的波函数。描绘电子运动的Schrodinger方程为:在式中，定义了一个能量E(R)，它是体系中固定原子核坐标时给定电子状态下的电子本征能量函数，称为分子势能函数或分子内势能函数(intramolecularpotentialenergyfunction)。在分子力场中，通常我们把E(R)称为位能面，或称Born-Oppenheimer面。假如一个解析表达式能拟合这个位能面，则此解析表达式就称为分子力场，简称力场。即力场是:描绘分

14、子构造和能量之间的一种数据解析表达式，包括原子核的坐标和一些可调节的参数。所以，一个力场确实定就是选择解析函数形式及确定参数2).力场能的表示形式分子力学的基本理论就是一个分子力场由分子内互相作用和分子间互相作用两大部分构成，分别对应于键伸缩、角弯曲、扭转运动、偶合互相作用以及Vdw互相作用和静电互相作用或氢键能。力场的参数和力常数可从一个分子迁移至类似构造的另一分子。一个力场中诸项的能量表达式基本上是类似的。其中Evalence，是各种成键互相作用能。成键互相作用能由键伸缩、角弯曲、二面角扭转和键角面外弯曲能等子项组成。即是:其中Ecrossterm是键、角之间的偶合(Coupling),所

15、相应的能量在力场中称为穿插项(CrossTerms)。这是由于考虑到分子内键长及键角的变化遭到相邻键、角的影响。包括键键偶合、键角偶合等。其中Enonbond是非键互相作用，由范德华能、静电效应及氢键能部分组成。详细到各种分子构造和分子势能函数的计算，由于研究的对象和考虑的因素不同而其解析表达式有所不同3).常用的分子力场介绍A.第一代经典的分子力场。主要是AMBER,CHARMm，及CVFF。其特点是:函数形式简单;应用的范围比拟特定(大部分合适于生物分子);优化力场参数的方法比拟多，结果也比拟好;能合理地预测分子构造，构象性质，凝聚态性质。B，第二代分子力场。CFF91,CFF,PCFF,

16、COMPASS及MMFF94。这一类力场的共同特点是函数形式较复杂，附加项多、适用范围宽，十分是大而复杂的分子模型的多能量极小和势垒、优化得到的力场参数比拟合理;能比拟好地预测分子构造，振动频率，构象性质。C.针对周期表中的所有元素的力场。所有的力学参数基于元素，杂化和化合的规则而产生，并为很多构造形式所证明。这类力场有:ESFF,UFF-VALBOND，DreidingFF。特点是:函数形式简单;适用的范围广;能比拟合理地预测分子构造;根照某些规则使力场参数化。Universal力场(UFF)是针对整个周期表的分子力学和动态模拟的力场，力场参数是依靠元素、元素的杂化及化合性而计算出来的，与电

17、荷平衡计算法相结合。D.用于特定目的的分子力场。如用于优化玻璃质的Glassforcefiel用于聚合体的PCFF，沸石构造的吸附作用的CVFF等。这类的特点自然是只局限于某一特定的范围内，不能在其他场合中使用。十.AMBER练习：真实水溶剂环境中的腺嘌呤胸腺嘧啶十聚体(A-DNA)的分子动力学模拟a)用nucgen命令创立A-DNA双螺旋构造i.建立一个a-dna.nucgen.in文件NUC1DA5AAAAAAAAA3NUC2DT5TTTTTTTTT3END$ABDNAb)运行nucgen命令创立nuc.pdb文件$AMBEROME/exe/nucgen-O-ia-dna.nucgen.i

18、n-oa-dna.nucgen.out-d$AMBERHOME/dat/leap/parm/nucgen.dat-pa-dna.nucgen.pdb这样产生a-dna.nucgen.out,a-dna.nucgen.pdb个文件。c)创立sander输入文件下一步是为了给sander创立拓扑文件(prmtop)和坐标文件(inpcrd)。在用xleap创立之前，我们需要给a-dna.nucgen.pdb文件中的第10、11个残基之间加TER标记。修改后的pdb文件为a-dna.pdb运行xleap：$AMBERHOME/exe/xleap-s-f$AMBERHOME/dat/leap/cmd/

19、leaprc.ff99在xleap窗口创立一个工作单元adna打开a-dna.pdb文件：adna=loadpdba-dna.pdb给我们的系统加抗衡粒子中性体系：addionsadnaNa+0如今产生中性体系(非溶剂化体系)的拓扑和坐标文件：saveamberparmadnaa-dna_charges_only.prmtopa-dna_charges_only.inpcrd下一步用一个截面八面体的水盒子使系统溶剂化：solvateoctadnaWATBOX2168.0#amber7.0不支持TIP3PBOX保存溶剂化后的拓扑和坐标文件：saveamberparmadnaa-dna_wat.p

20、rmtopa-dna_wat.inpcrdd)对系统进行能量优化：在对系统做MD之前要对系统做一次能量优化，为了除去加氢和溶剂化时传生的不良触点。第一步：优化离子和水的位置。优化的前500步用最陡下降法后500步用共轭梯度法，在DNA分子(第120个残基)上采用500kal/mol的限制作用力。a-dna_min1.inA-DNA10-mer:initialminimisationsolvent+ions&cntrlimin=1,maxcyc=1000,ncyc=500,ntb=1,ntr=1,cut=10/HoldtheDNAfixed500.0RES120ENDENDsander-O-ia

21、-dna_min1.in-oa-dna_min1.out-pa-dna_wat.prmtop-ca-dna_wat.inpcrd-ra-dna_min1.rst-refa-dna_wat.inpcrd第二步：全局优化。前1000步采用最陡下降法，后1500步采用共轭梯度法，这时不需要约束力。输入文件如下：a-dna_min2.inA-DNA10-mer:initialminimisationsolvent+ions&cntrlimin=1,maxcyc=2500,ncyc=1000,ntb=1,ntr=0,cut=10/sander-O-ia-dna_min2.in-oa-dna_min2.o

22、ut-pa-dna_wat.prmtop-ca-dna_min1.rst-ra-dna_min2.rste)初始平衡这一步开场运行分子动力学。我们首先运行20ps的恒V加热平衡。采用的方法是：在DNA上加10kal/mol的弱限制力，采用恒V周期边界条件做20ps的MD，起始温度是0K，最终温度是300K，采用Langevin温度控制，在H原子上采用SHAKE限制，步长为2fs。adna_md1.inA-DNA10-mer:20psMDwithresonDNA&cntrlimin=0,irest=0,ntx=1,ntb=1,cut=10,ntr=1,ntc=2,ntf=2,tempi=0.0,

23、temp0=300.0,ntt=3,gamma_ln=1.0,nstlim=10000,dt=0.002,ntpr=250,ntwx=250,ntwr=10000/KeepDNAfixedwithweakrestraints10.0RES120ENDENDsander-O-ia-dna_md1.in-oa-dna_md1.out-pa-dna_wat.prmtop-ca-dna_min2.rst-ra-dna_md1.rst-refa-dna_min2.rst-xa-dna_md1.mdcrdf)进一步的平衡这一步将运行更长的1.8ns的MD，对DNA不加以限制，采用Langevin恒温法使温

24、度保持在300K，采用恒压周期边界条件(1atm)，对H原子采用SHAKE限制，步长为2fs。adna_md_1800ps.inA-DNA10-mer:1800psofMD&cntrlimin=0,irest=1,ntx=7,ntb=2,pres0=1.0,ntp=1,taup=2.0,cut=10,ntr=0,ntc=2,ntf=2,tempi=300.0,temp0=300.0,ntt=3,gamma_ln=1.0,nstlim=900000,dt=0.002,ntpr=250,ntwx=250,ntwr=10000/$AMBERHOME/exe/sander-O-ia-dna_md_18

25、00ps.in-oa-dna_md2.out-pa-dna_wat.prmtop-ca-dna_md1.rst-ra-dna_md2.rst-xa-dna_md2.mdcrdg)结果分析：从输出文件a-dna_md1.outa-dna_md2.out中能够看书系统的很多特性，如：系统的势能、动能、总能量、温度、压力、体积、密度、均方根差。先用process-mdout.perl脚本对a-dna_md2.out进行数据处理，产生一系列summary文件，再用绘图程序xmgrace-5.99.0对上面产生的summary文件绘图。1.通过对势能、动能、总能量、温度、压力、体积、密度的分析，在系统得

26、到最初的释放后，势能、动能、总能量、温度、压力、体积、密度保持恒定。这讲明在这次模拟中没有出现大的问题。2.接下来通过计算主链上原子的RMSd来看看DNA结果发生了什么变化。使用ptraj来对整个分子轨道重新成像并除去分子轨道文件里所有的水。这样做的原因是不是要在VMD显示中隐藏水，而是如今的mdcrd文件要比先前的大得多，假如把所有的东西都参加到VMD中会使VMD运行很慢。在分子动力学模拟中，我们只对DNA主链感兴趣，所以只关注主链原子P,O3,O5,C3,C4,C5。combine_mdcrds_and_strip.ptrajtrajina-dna_md1.mdcrdtrajina-dna_md2.mdcrdtrajouta-dna_0-1820ps_no_wat.mdcrdrmsfirstouta-dna_0-1820ps.rmsfitP,O3,O5,C3,C4,C5time0.50center:1-20imagefamiliarstrip:WATptraja-dna_wat.prmtop