高分子材料纳米级粘着接触及分子动力学方法-讲座资料2010课件

高分子材料纳米级粘着接触及分子动力学模拟方法计算材料科学中的模型、算法和多尺度关联系列讲座与讨论西北工业大学科研处西北工业大学理学院西北工业大学材料学院联合举办西北工业大学高性能计算中心高分子材料纳米级粘着接触及III、纳米、亚纳米尺度一、第一原理方法及其在材料科学中的应用二、密度泛函方法及其应用三、从原子分子到纳米尺度

——分子力学、分子动力学方法及其应用III、纳米、亚纳米尺度从原子分子到纳米尺度——分子力学、分子动力学方法及其应用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年代以来得到迅猛发展和广泛应用一、MM、MD基础理论基本思想事先构造出简单体系(如链段、官能团等各种不同结构的小片段)的势能函数,简称势函数或力场(forcefield)将势函数建成数据库，在形成较大分子的势函数时，从数据库中检索到结构相同的片段，组合成大体系的势函数利用分子势能随原子位置的变化有极小值的性质，确定大分子的结构即为分子力学(MM)利用势函数，建立并求解与温度和时间有关的牛顿运动方程，得到一定条件下体系的结构随时间的演化关系即为分子动力学(MD)理论方法的核心是构造势函数势函数：势能与原子位置的关系。且往往是不知道的需要通过其他方法，如量子化学方法及实验数据获得基本思想rrθθrE过渡状态

E反应物产物分子势函数曲面——势能面示意图rrθθrE过渡状态E反应物产物分子势函数曲面——势能面示1、分子力场分子片段力场的函数表达式中包含自变量和力场参数其中自变量为分子的结构参数，独立参数为键长、键角和二面角，如图而力场参数一般通过与实验数据和从头算数据进行最小二乘法拟合来确定b’bθbθθ’有的还使用一个非独立参数：面外弯曲角χφ1、分子力场而力场参数一般通过与实验数据和从头算数势函数形式很多，目前已被广泛使用的力场有如CFF、MM2、MM3、MM4、MMFF、AMBER、CHARMM、DREIDING、UFF和COMPASS等形式虽多,但一般总表达为分子内与分子间势能之和:V总=V键合+V非键合分子内势能(键合)包括键伸缩、键角弯曲和二面角扭转势能分子间势能(非键合)包括范德华势和静电势,有的还包括H键：V键合=V键伸缩+V键角弯曲+V二面角扭转V非键合=V范德华+V静电+V氢键键合势函数中，一些力场还包含交叉项，使精度更高交叉项的含义：如键长变化时，键角弯曲势能随键长的不同而不同，等势函数形式很多，目前已被广泛使用的力场有如CFF、MM2、M例:

COMPASS-98力场(condensed-phaseoptimizedmolecularpotentialsforatomisticsimulationstudies)的表达式如下每个k是一独立的力场参数,下标“0”代表参考(~平衡)结构参数：显然是一个非谐性力场例:COMPASS-98力场(condensed-phas力场参数{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}力场参数{k}最小二乘法确定非键合势函数中，静电相互作用表示分子中各原子静电荷的库仑相互作用对势能的贡献不同的力场，静电相互作用表达式基本相同范德华势也大都采用Lennard-Jones函数，但函数中的指数有所不同。如COMPASS-98的非键合势函数为：静电相互作用：范德华势：即Lennard-Jones9-6函数其它力场范德华势较多采用L-J12-6函数：非键合势函数中，静电相互作用表示分子中各原子静电荷的库仑相互当然，在建立分子的势能函数时，还有一些更细致的问题要考虑，如：势能展开项的截断周期结构的处理多组分混合物体系含有离子的体系或金属中的离子等不同的方法或程序中，分别都有更详细的讨论。参考：[德]D.罗伯.计算材料学.北京：化学工业出版社,2002,9俞庆森,朱龙观.分子设计导论.北京：高等教育出版社,2000杨小震.分子模拟与高分子材料.北京：科学出版社，2002熊家炯主编.材料设计.天津：天津大学出版社，2000SunH，RenP,FriedJR.TheCOMPASSForceField:ParameterizationandValidationforpolyphosphazenes.ComputationalandTheoreticalPolymerScience,1998,8(1/2):229SunH.COMPASS:AnabInitioForce-FieldOptimizedforCondensed-PhaseApplications-OverviewwithDetailsonAlkaneandBenzeneCompounds.J.Phys.Chem.,1998,102:7338

当然，在建立分子的势能函数时，还有一些更细致的问题要考虑，如2、分子力学方法MM是确定分子结构的方法利用分子势能随结构的变化而变化的性质，确定分子势能极小时的平衡结构(stationarypoint)物理模型:视原子为质点，视化学键为弹簧，而弹力常数完全由数据库中的分子力场来确定因此是直接用势函数研究问题，不考虑原子的动能不考虑动能所对应的结构，相当于体系处于T=0K时的结果2、分子力学方法由力场首先构造并得到分子的势函数E(x)利用在E(x)的极小点处，E(x)随各原子独立的空间坐标{xi,i=1,2,3,…3N-6}(与内坐标{ri、θi、φi}等价)的一阶微分=0以及全部二阶微分>0的数学条件：(i=1,2,3,…3N-6)进行结构优化具体步骤：

1)用各种方法构造出一个分子的任意结构，得到初始的结构参数{xi0,i=1,2,3,…3N-6}2)进行坐标变换，即根据需要将原子的直角坐标转变成内坐标{ri、θi、φi}或反之由力场首先构造并得到分子的势函数E(x)(i=1,2,3,…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计算中，有时还要考虑体系所处的外压条件，此时，可通过压力因子的定义来调节原子坐标(i=1,2,3,…3N-6)1)用各种方法构造出一个分子的任意结构，得到初始的结构参数3、分子动力学方法分子力场是分子的静态势函数而实际过程通常是在一定温度和一定压力下发生的为了更切实际地了解体系运动和演化的过程，必须考虑体系中原子的运动，并与温度T和时间t建立联系我们知道，温度是原子分子热运动剧烈程度的量度根据统计热力学，对于n个原子的体系，体系的温度T与各原子的运动速率vi的关系为：又因体系中各原子的速率为vi时，动量pi=mivi，对应总动能K(p)为：3、分子动力学方法又因体系中各原子的速率为vi时，动量p势能由力场确定为E(x)，因此体系的Hamilton量H为：与量子力学不同，经典力学对Hamilton量不进行算符化处理，也不建立和求解本征方程，而是建立并求解经典运动方程：计算过程一般为：在一定的统计系综下1)由原子位置和连接方式，从数据库调用力场参数并形成体系势函数2)由给定温度计算体系动能以及总能量3)计算各原子的势能梯度,得到原子在力场中所受的力即dp/dt=mdv/dt=ma=F势能由力场确定为E(x)，因此体系的Hamilton量H为：1)由原子位置和连接方式，调用力场参数并形成体系势函数2)由给定温度计算体系动能以及总能量3)计算各原子的势能梯度,得到原子在力场中所受的力

即dp/dt=mdv/dt=ma=F4)对每个原子，在一定时间间隔内，用牛顿方程求解其运动行为：5)显示体系能量和结构6)取下一时间间隔，返回步骤1)不断循环反复，可设定循环次数或强行终止计算其中,不仅时间间隔可以根据需要取不同大小(一般~1fs=10-15s),温度可以任意设定，而且还可以在循环过程中逐渐改变温度，即研究体系的退火(annealing)行为1)由原子位置和连接方式，调用力场参数并形成体系势函数有时需要进一步考虑外场的作用，如压力、电场、磁场、重力场等从原则上讲，这些问题都不难解决，而且还在进一步发展中，不再赘述此外，即使是分子力学与分子动力学方法，也受计算量的限制，所处理的体系不可能太大，好的计算机可达到数十万~百万个原子的规模(20~50nm以下/MD106~108个原子——罗伯p110)对于更大(>50nm、微米)尺度的问题，人们也在发展介观尺度(mesoscale)的方法，并且已取得一些成果分子动力学方法也有一些本质的缺陷如：势函数精度(客观性)势函数形式在每次计算中都不变,故不能模拟如分子在高温下,结构发生断裂的热裂解过程(国内引入开关函数已开始做)此时,可用QM处理化学变化或e转移,余用MM，称QM/MM有时需要进一步考虑外场的作用，如压力、电场、磁场、重力场等二、MM、MD计算程序因程序结构相对简单，国内外有不少,且还在不断产生如：Chem3D,ChemOffice,Cerius2,

MaterialStudio,Alchemy2000,Sybyl,Biosym,HyperChemx.x,Spartanx.x,Chemgraf,Bilder,Script,COGS,Gaussian-03;MP(MolecularProperties,杨小震)等输入：体系模型建模(Builder)，选择力场、系综MM:无特殊输入MD:循环次数(万~几10万)、温度、时间间隔(fs),外场,温度变化速率，其它性质等输出：MM:势能-结构曲线(数值)，动态结构图,其它性质MD:(同上)，及其它MD性质注：第一原理及MM一般只能得到势能极小点结构，而MD可越过一些小的势垒，甚至可达到最小点二、MM、MD计算程序三、MM、MD方法的应用领域：高分子、生命科学、药物设计、催化、半导体其它功能材料、结构材料等分子力学是用计算机在原子水平上模拟给定分子模型的结构与性质，进而得到分子的各种物理性质与化学性质，如结构参数、振动频率、构象能量、相互作用能量、偶极矩、密度、摩尔体积、汽化焓等分子动力学方法能实时将分子的动态行为显示到计算机屏幕上,便于直观了解体系在一定条件下的演变过程MD含温度与时间,因此还可得到如材料的玻璃化转变温度、热容、晶体结晶过程、输送过程、膨胀过程、动态弛豫(relax)以及体系在外场作用下的变化过程等三、MM、MD方法的应用分子模拟较早应用于高分子问题的研究应用范围主要包括:物理性质、结构、构象与弹性、晶体结构、力学性能、玻璃态与玻璃化转变、光谱性质、非线性光学性质、电性质、共混与分子间相互作用等应用于生物科学和药物设计也十分普及，如蛋白质的多级结构与性质，病毒、药物作用机理、特效药物的大通量筛选与快速开发等在化学领域，用于表面催化与催化机理、溶剂效应、原子簇的结构与性质研究等在材料科学,用于材料的优化设计、结构与力学性能、热加工性能预报、界面相互作用、纳米材料结构与性能研究等分子模拟较早应用于高分子问题的研究罗伯，计算材料学，p129:已有很多研究者采用分子动力学在原子尺度上进行了材料性质及结构的模拟在材料科学领域，相变及晶格缺陷结构的研究一直是普遍关注的热点在这方面，人们关注的焦点集中在同相和异相界面、位错、裂纹、界面偏析、失配(misfit)位错、次级晶界位错以及聚合物中的结构筹等问题而且，原子尺度上的材料合成和设计已经取得实质性进展也正是在这一点上，应力被分别地引入到金属、聚合物和陶瓷材料的具体研究当中已发表的工作包括：经典例子和基本原理、外界环境条件、失配位错的结构、界面结构及其迁移率、晶格位错的结构和能量、偏析效应、微裂纹结构、薄膜和表面、异相界面、材料设计、原子论和有限元罗伯，计算材料学，p129:谢谢大家！谢谢大家！高分子材料纳米级粘着接触及分子动力学模拟方法计算材料科学中的模型、算法和多尺度关联系列讲座与讨论西北工业大学科研处西北工业大学理学院西北工业大学材料学院联合举办西北工业大学高性能计算中心高分子材料纳米级粘着接触及III、纳米、亚纳米尺度一、第一原理方法及其在材料科学中的应用二、密度泛函方法及其应用三、从原子分子到纳米尺度

——分子力学、分子动力学方法及其应用III、纳米、亚纳米尺度从原子分子到纳米尺度——分子力学、分子动力学方法及其应用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年代以来得到迅猛发展和广泛应用一、MM、MD基础理论基本思想事先构造出简单体系(如链段、官能团等各种不同结构的小片段)的势能函数,简称势函数或力场(forcefield)将势函数建成数据库，在形成较大分子的势函数时，从数据库中检索到结构相同的片段，组合成大体系的势函数利用分子势能随原子位置的变化有极小值的性质，确定大分子的结构即为分子力学(MM)利用势函数，建立并求解与温度和时间有关的牛顿运动方程，得到一定条件下体系的结构随时间的演化关系即为分子动力学(MD)理论方法的核心是构造势函数势函数：势能与原子位置的关系。且往往是不知道的需要通过其他方法，如量子化学方法及实验数据获得基本思想rrθθrE过渡状态

E反应物产物分子势函数曲面——势能面示意图rrθθrE过渡状态E反应物产物分子势函数曲面——势能面示1、分子力场分子片段力场的函数表达式中包含自变量和力场参数其中自变量为分子的结构参数，独立参数为键长、键角和二面角，如图而力场参数一般通过与实验数据和从头算数据进行最小二乘法拟合来确定b’bθbθθ’有的还使用一个非独立参数：面外弯曲角χφ1、分子力场而力场参数一般通过与实验数据和从头算数势函数形式很多，目前已被广泛使用的力场有如CFF、MM2、MM3、MM4、MMFF、AMBER、CHARMM、DREIDING、UFF和COMPASS等形式虽多,但一般总表达为分子内与分子间势能之和:V总=V键合+V非键合分子内势能(键合)包括键伸缩、键角弯曲和二面角扭转势能分子间势能(非键合)包括范德华势和静电势,有的还包括H键：V键合=V键伸缩+V键角弯曲+V二面角扭转V非键合=V范德华+V静电+V氢键键合势函数中，一些力场还包含交叉项，使精度更高交叉项的含义：如键长变化时，键角弯曲势能随键长的不同而不同，等势函数形式很多，目前已被广泛使用的力场有如CFF、MM2、M例:

COMPASS-98力场(condensed-phaseoptimizedmolecularpotentialsforatomisticsimulationstudies)的表达式如下每个k是一独立的力场参数,下标“0”代表参考(~平衡)结构参数：显然是一个非谐性力场例:COMPASS-98力场(condensed-phas力场参数{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}力场参数{k}最小二乘法确定非键合势函数中，静电相互作用表示分子中各原子静电荷的库仑相互作用对势能的贡献不同的力场，静电相互作用表达式基本相同范德华势也大都采用Lennard-Jones函数，但函数中的指数有所不同。如COMPASS-98的非键合势函数为：静电相互作用：范德华势：即Lennard-Jones9-6函数其它力场范德华势较多采用L-J12-6函数：非键合势函数中，静电相互作用表示分子中各原子静电荷的库仑相互当然，在建立分子的势能函数时，还有一些更细致的问题要考虑，如：势能展开项的截断周期结构的处理多组分混合物体系含有离子的体系或金属中的离子等不同的方法或程序中，分别都有更详细的讨论。参考：[德]D.罗伯.计算材料学.北京：化学工业出版社,2002,9俞庆森,朱龙观.分子设计导论.北京：高等教育出版社,2000杨小震.分子模拟与高分子材料.北京：科学出版社，2002熊家炯主编.材料设计.天津：天津大学出版社，2000SunH，RenP,FriedJR.TheCOMPASSForceField:ParameterizationandValidationforpolyphosphazenes.ComputationalandTheoreticalPolymerScience,1998,8(1/2):229SunH.COMPASS:AnabInitioForce-FieldOptimizedforCondensed-PhaseApplications-OverviewwithDetailsonAlkaneandBenzeneCompounds.J.Phys.Chem.,1998,102:7338

当然，在建立分子的势能函数时，还有一些更细致的问题要考虑，如2、分子力学方法MM是确定分子结构的方法利用分子势能随结构的变化而变化的性质，确定分子势能极小时的平衡结构(stationarypoint)物理模型:视原子为质点，视化学键为弹簧，而弹力常数完全由数据库中的分子力场来确定因此是直接用势函数研究问题，不考虑原子的动能不考虑动能所对应的结构，相当于体系处于T=0K时的结果2、分子力学方法由力场首先构造并得到分子的势函数E(x)利用在E(x)的极小点处，E(x)随各原子独立的空间坐标{xi,i=1,2,3,…3N-6}(与内坐标{ri、θi、φi}等价)的一阶微分=0以及全部二阶微分>0的数学条件：(i=1,2,3,…3N-6)进行结构优化具体步骤：

1)用各种方法构造出一个分子的任意结构，得到初始的结构参数{xi0,i=1,2,3,…3N-6}2)进行坐标变换，即根据需要将原子的直角坐标转变成内坐标{ri、θi、φi}或反之由力场首先构造并得到分子的势函数E(x)(i=1,2,3,…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计算中，有时还要考虑体系所处的外压条件，此时，可通过压力因子的定义来调节原子坐标(i=1,2,3,…3N-6)1)用各种方法构造出一个分子的任意结构，得到初始的结构参数3、分子动力学方法分子力场是分子的静态势函数而实际过程通常是在一定温度和一定压力下发生的为了更切实际地了解体系运动和演化的过程，必须考虑体系中原子的运动，并与温度T和时间t建立联系我们知道，温度是原子分子热运动剧烈程度的量度根据统计热力学，对于n个原子的体系，体系的温度T与各原子的运动速率vi的关系为：又因体系中各原子的速率为vi时，动量pi=mivi，对应总动能K(p)为：3、分子动力学方法又因体系中各原子的速率为vi时，动量p势能由力场确定为E(x)，因此体系的Hamilton量H为：与量子力学不同，经典力学对Hamilton量不进行算符化处理，也不建立和求解本征方程，而是建立并求解经典运动方程：计算过程一般为：在一定的统计系综下1)由原子位置和连接方式，从数据库调用力场参数并形成体系势函数2)由给定温度计算体系动能以及总能量3)计算各原子的势能梯度,得到原子在力场中所受的力即dp/dt=mdv/dt=ma=F势能由力场确定为E(x)，因此体系的Hamilton量H为：1)由原子位置和连接方式，调用力场参数并形成体系势函数2)由给定温度计算体系动能以及总能量3)计算各原子的势能梯度,得到原子在力场中所受的力

即dp/dt=mdv/dt=ma=F4)对每个原子，在一定时间间隔内，用牛顿方程求解其运动行为：5)显示体系能量和结构6)取下一时间间隔，返回步骤1)不断循环反复，可设定循环次数或强行终止计算其中,不仅时间间隔可以根据需要取不同大小(一般~1fs=10-15s),温度可以任意设定，而且还可以在循环过程中逐渐改变温度，即研究体系的退火(annealing)行为1)由原子位置和连接方式，调用力场参数并形成体系势函数有时需要进一步考虑外场的作用，如压力、电场、磁场、重力场等从原则上讲，这些问题都不难解决，而且还在进一步发展中，不再赘述此外，即使是分子力学与分子动力学方法，也受计算量的限制，所处理的体系不可能太大，好的计算机可达到数十万~百万个原子的规模(20~50nm以下/MD106~108个原子——罗伯p110)对于更大(>50nm、微米)尺度的问题，人们也在发展介观尺度(mesoscale)的方法，并且已取得一些成果分子动力学方法也有一些本质的缺陷如：势函数精度(客观性)势函数形式在每次计算中都不变,故不能模拟如分子在高温下,结构发生断裂的热裂解过程(国内引入开关函数已开始做)此时,可用QM处理化学变化或e转移,余用MM，称QM/MM有时需要进一步考虑外场的作用，如压力、电场、磁场、重力场等二、MM、