分子力学和分子动力学

学习目标了解计算化学的主要方法和基本概念能够读懂相关领域的文献，了解常见分子模拟软件及使用掌握基本的计算化学研究思路和方法。了解分子模拟能干什么，我们可以作那些工作，作到什么精度，为在自己的研究领域开展分子模拟方面的研究打一个基础。分子模拟的双重性质分子模拟具有理论和实验的双重性质分子模拟不能完全取代实验理论实验模拟理论的正确性模拟参数的正确性模拟方法的选择理论的更新系综•系综（ensemble）代表一大群相类似的体系的集合。对一类相同性质的体系，其微观状态（比如每个粒子的位置和速度）仍然可以大不相同。（实际上，对于一个宏观体系，所有可能的微观状态数是天文数字。）统计物理的一个基本假设（各态历经假设）是：对于一个处于平衡的体系，物理量的时间平均，等于对对应系综里所有体系进行平均的结果。体系的平衡态的物理性质可以对不同的微观状态求和来得到。系综的概念是由约西亚·威拉德·吉布斯（J.WillardGibbs）在1878年提出的。常用系综微正则系综(microcanonicalensemble)正则系综(canonicalensemble)巨正则系综(grandcanonicalensemble)等温等压系综(isothermal-isobaricensemble)常用系综•微正则系综(microcanonicalensemble)：系综里的每个体系具有相同的能量（通常每个体系的粒子数和体积也是相同的）。•正则系综(canonicalensemble)：系综里的每个体系都可以和其他体系交换能量（每个体系的粒子数和体积仍然是固定且相同的），但是系综里所有体系的能量总和是固定的。系综内各体系有相同的温度。•巨正则系综(grandcanonicalensemble)：正则系综的推广，每个体系都可以和其他体系交换能量和粒子，但系综内各体系的能量总和以及粒子数总和都是固定的。（系综内各体系的体积相同。）系综内各个体系有相同的温度和化学势。•等温等压系综(isothermal-isobaricensemble)：正则系综的推广，体系间可交换能量和体积，但能量总和以及体积总和都是固定的。（系综内各体系有相同的粒子数。）正如它的名字，系综内各个体系有相同的温度和压强。分子动力学方法工作框图给定t时刻的坐标和速度以及其他动力学信息，那么就可计算出t＋Δt时刻的坐标和速度。程序构成方式•①输入指定运算条件的参数（初始温度，粒子数，密度，时间步长）；•②体系初始化（选定初始坐标和初始速度）；•③计算作用在所有粒子上的力；•④解牛顿运动方程（第③和第④步构成了模拟的核心，重复这两步，直到体系的演化到指•定的时间）；•⑤计算并输出物理量的平均值，完成模拟。积分算法优劣的判据分子动力学中一个好的积分算法的判据主要包括：•①计算速度快；•②需要较小的计算机内存；•③允许使用较长的时间步长；•④表现出较好的能量守恒。分子动力学的适用范围•分子动力学方法只考虑多体系统中原子核的运动，而电子的运动不予考虑，量子效应忽略。经典近似在很宽的材料体系都较精确；但对于涉及电荷重新分布的化学反应、键的形成与•断裂、解离、极化以及金属离子的化学键都不适用，此时需要使用量子力学方法。经典分子动•力学方法（MD）也不适用于低温，因为量子物理给出的离散能级之间的能隙比体系的热能大，•体系被限制在一个或几个低能态中。当温度升高或与运动相关的频率降低（有较长的时间标•度）时，离散的能级描述变得不重要，在这样的条件下，更高的能级可以用热激发描述。•电子运动具有更高的特征频率，必须用量子力学以及量子／经典理论联合处理。这些技术近年来取得了很大进步。在这些方法中，体系中化学反应部分用量子理论处理，而其他部分用经典模型处理。QM/MM分子力学、分子动力学方法及其应用MolecularMechanics/MolecularDynamics•一、MM、MD理论基础•二、MM、MD计算程序•三、MM、MD方法的应用一、MM、MD基础理论•原则上,第一原理方法在理论上已经能解决所有问题•但计算量太大，计算机资源有限，原子数目较多时，如高分子、蛋白质、原子簇以及研究表面问题、功能材料或材料的力学性能等，实际上难以完成计算•为此，发展了分子力学(MolecularMechanics,MM)与分子动力学(MolecularDynamics,MD)方法•它们的应用，又称分子模拟(molecularsimulation,molecularmodeling)或分子设计(moleculardesign)•MM与MD是经典力学方法，针对的最小结构单元不再是电子而是原子•因原子的质量比电子大很多，量子效应不明显，可近似用经典力学方法处理•20世纪30年代,Andrews最早提出分子力学(MM)的基本思想；40年代以后得到发展,并用于有机小分子研究。90年代以来得到迅猛发展和广泛应用基本思想•事先构造出简单体系(如链段、官能团等各种不同结构的小片段)的势能函数,简称势函数或力场(forcefield)•将势函数建成数据库，在形成较大分子的势函数时，从数据库中检索到结构相同的片段，组合成大体系的势函数•利用分子势能随原子位置的变化有极小值的性质，确定大分子的结构即为分子力学(MM)•利用势函数，建立并求解与温度和时间有关的牛顿运动方程，得到一定条件下体系的结构随时间的演化关系即为分子动力学(MD)•理论方法的核心是构造势函数•势函数：势能与原子位置的关系。且往往是不知道的•需要通过其他方法，如量子化学方法及实验数据获得rrθθrE过渡状态E反应物产物分子势函数曲面——势能面示意图•1、分子力场•分子片段力场的函数表达式中包含自变量和力场参数•其中自变量为分子的结构参数，独立参数为键长、键角和二面角，如图•而力场参数一般通过与实验数据和从头算数据进行最小二乘法拟合来确定b’bθbθθ’有的还使用一个非独立参数：面外弯曲角χφ•势函数形式很多，目前已被广泛使用的力场有如CFF、MM2、MM3、MM4、MMFF、AMBER、CHARMM、DREIDING、UFF和COMPASS等•形式虽多,但一般总表达为分子内与分子间势能之和:V总=V键合+V非键合•分子内势能(键合)包括键伸缩、键角弯曲和二面角扭转势能•分子间势能(非键合)包括范德华势和静电势,有的还包括H键：V键合=V键伸缩+V键角弯曲+V二面角扭转V非键合=V范德华+V静电+V氢键•键合势函数中，一些力场还包含交叉项，使精度更高•交叉项的含义：如键长变化时，键角弯曲势能随键长的不同而不同，等•例:COMPASS-98力场(condensed-phaseoptimizedmolecularpotentialsforatomisticsimulationstudies)的表达式如下•每个k是一独立的力场参数,下标“0”代表参考(~平衡)结构参数：bbbbkbbkbbkE])()()([404303202键伸缩：])()()([404303202kkkE键弯曲：)]}3cos(1[)2cos(1[)]cos(1[{3,032,021,01kkkE二面角：2kE键面外弯曲：',,'0'0'0'0))(())((bbbbxkbbbbkE交叉项：,,321000]3cos2coscos)[())((bbbkkkbbbbk,'321'0']3cos2coscos)[(bkkkbb,3210]3cos2coscos)[(kkk',,'0'0'))((cosk•力场参数{k}最小二乘法确定基本思想•如：R-COOH基团•1)由abinitio(构型优化方法)计算出平衡结构，得到参考结构参数{bi0,θi0,φi0}•2)用伪随机数方法将{bi0,θi0,φi0}人为改变成若干(n)组非平衡结构参数{bi,θi,φi}l(l=1,2,…,n)•3)用各{bi,θi,φi}l结构参数分别进行abinitio计算，得到{bi,θi,φi}l结构参数下对应的能量El(l=1,2,…,n)•4)将El(l=1,2,…,n)和{bi,θi,φi}l代入势能表达式•5)用最小二乘法拟合，确定力场参数{k}COOH•非键合势函数中，静电相互作用表示分子中各原子静电荷的库仑相互作用对势能的贡献•不同的力场，静电相互作用表达式基本相同•范德华势也大都采用Lennard-Jones函数，但函数中的指数有所不同。如COMPASS-98的非键合势函数为：静电相互作用：jiijjirqqEelec范德华势：jiijijijijijvdWrrrrE,609032即Lennard-Jones9-6函数其它力场范德华势较多采用L-J12-6函数：jiijijijijijvdWrrrrDE,6012002)(•当然，在建立分子的势能函数时，还有一些更细致的问题要考虑，如：–势能展开项的截断–周期结构的处理–多组分混合物体系–含有离子的体系或–金属中的离子等•不同的方法或程序中，分别都有更详细的讨论。参考：•[德]D.罗伯.计算材料学.北京：化学工业出版社,2002,9•俞庆森,朱龙观.分子设计导论.北京：高等教育出版社,2000•杨小震.分子模拟与高分子材料.北京：科学出版社，2002•熊家炯主编.材料设计.天津：天津大学出版社，2000•SunH，RenP,FriedJR.TheCOMPASSForceField:ParameterizationandValidationforpolyphosphazenes.ComputationalandTheoreticalPolymerScience,1998,8(1/2):229•SunH.COMPASS:AnabInitioForce-FieldOptimizedforCondensed-PhaseApplications-OverviewwithDetailsonAlkaneandBenzeneCompounds.J.Phys.Chem.,1998,102:7338•2、分子力学方法•MM是确定分子结构的方法•利用分子势能随结构的变化而变化的性质，确定分子势能极小时的平衡结构(stationarypoint)•物理模型:视原子为质点，视化学键为弹簧，而弹力常数完全由数据库中的分子力场来确定•因此是直接用势函数研究问题，不考虑原子的动能•不考虑动能所对应的结构，相当于体系处于T=0K时的结果•由力场首先构造并得到分子的势函数E(x)•利用在E(x)的极小点处，E(x)随各原子独立的空间坐标{xi,i=1,2,3,…3N-6}(与内坐标{ri、θi、φi}等价)的一阶微分=0以及全部二阶微分0的数学条件：0)(0)(22iixxExxE(i=1,2,3,…3N-6)进行结构优化具体步骤：•1)用各种方法构造出一个分子的任意结构，得到初始的结构参数{xi0,i=1,2,3,…3N-6}•2)进行坐标变换，即根据需要将原子的直角坐标转变成内坐标{ri、θi、φi}或反之•1)用各种方法构造出一个分子的任意结构，得到初始的结构参数{xi0,i=1,2,3,…3N-6}•2)进行坐标变换，即根据需要将原子的直角坐标转变成内坐标(键长、键角、二面角)或反之•3)建立分子体系的势能表达式E(x)•4)计算E(x)随各坐标的一阶、二阶导数•5)计算接近数学条件的坐标增量•6)得到新的结构参数{xi1,i=1,2,3,…3N-6}•重复4)、5)、6)，直至最后两次得到的体系势能之差或总体势能梯度的均方根值达到预定精度范围为止•MM计算中，有时还要考虑体系所处的外压条件，此时，可通过压力因子的定义来调节原子坐