1、第 1 章 前 言1.1 研究背景1.1.1 分子模拟及其发展分子模拟(Molecular Simulation)为二十世纪初发展起来的一种计算机模拟方法，它泛指用于模拟分子或分子体系性质的方法，主要用于探索研究具有三维结构的分子结构和分子的性能 i。分子模拟是根据物理和化学的基本原理构建一个模型（通常是数学模型,是对某种分子体系或反应过程的理想化描述) ，建立一种以计算数据(由计算机来执行)来代替实验测量的研究方法,并获取相关的物理和化学信息。分子模拟在材料科学方面的应用包括模拟材料的结构、计算材料的性质、预测材料的行为、验证实验结果（重现实验过程）、从微观角度认识材料，总之是为了更深层次理
2、解材料的结构，认识材料的各种行为。分子模拟的主要优势在于可以降低实验成本、具有较高的安全性、实现通常条件下较难或无法进行的实验（例如：超低温，低于-100；超高压，高于 100Mpa）、研究极快速的反应和变化等。 R.S.Mulliken 获得诺贝尔获奖时的感言是：“总之，我愿意强调我的信念：计算化学的年代已经到来，成千上百的化学家以计算机代替实验室，来获得众多的化学信息。唯一的障碍是你必须偿付机时费。” 阿基米德曾说过：“给我一个支点，我就能够翘起地球。”，但是分子模拟告诉我们:“Give me an enough powerful computer ,I can simulate the
3、whole world”。从 1980 年开始，每年在Engineering Village 中关于“Molecular Simulation”的文章数目由 37 篇递增到最高5209（2008 年）篇。与分子模拟有关的论文，美国（United States）发表的篇数最多，高达 16351 篇，其次是日本，中国名列第三。分子模拟作为一种计算机模拟技术，主要可以进行解释工作和预测工作。前者为实验奠定理论基础，通过模拟解释实验现象、建立理论模型、探讨过程机理等，后者为实验过程提供可能性和可行性研究，进行方案辅助设计、材料性能预测、过程优化筛选等。不同的分子模拟方法可以得到不同的信息。量子力学模拟
4、方法可以计算得到分子的大多数性质，如结构、构象、偶极矩、电离能、电子亲和力、电子密度、过渡态和反应途径等；分子力学可以计算分子体系的稳定构象、热力学特性、振动光谱等；能量最小化可以探索相空间（phase space)和势能面（ potential curve)，可以找出局部(local) 与全局(global)的最小点及转化过程的马鞍点(saddle point)；Monte Carlo 可以计算复杂分子体系的结构变化，特别是相变化；分子动力学可以得到复杂分子的热力学性质、结构、力学性质，特别是可以观察体系的动态演变，得到许多与时间有关的热动力学性质；布朗动力学可以研究蛋白质在水溶液中的折叠过
5、程；构象分析可以研究复杂分子稳态和亚稳态结构之间的演变等等。分子模拟所涉及的研究领域，涵盖了物理、化学、化工、材料、生化等几乎一切可以通过建立理论模型进行研究的体系，多数能够得到与实验结果近似的计算结果，所以分子模拟己经逐渐成为与实验技术并重的强有力的研究手段。分子模拟实际上并不仅仅局限于计算机模拟，但今天的分子模拟己和计算机模拟密不可分，正是由于计算机的高速计算技术的发展才使得分子模拟能够像今天这样发挥如此重要的作用。利用计算先行以了解更多的分子特性，已成为合成化学家和药物设计学家所依赖的重要方法。化学家们通过这种方法可以设计出最佳的反应途径，预测合成的可能性，且节省许多时间和避免材料的浪费
6、。分子模拟除了在药物设计方面的应用外，已被广泛地应用于研究金属材料、无机材料、高分子聚合物材料、生物材料等复杂庞大的体系。由于模拟结果随着计算机的发展与模拟方法（算法）的改良而更加精确，分子模拟已成为如今化学各个领域中不可缺少的工具。1.1.2 分子模拟的主要方法分子模拟的主要模拟方法有量子力学模拟和经典力学模拟，量子力学模拟主要根据从头算方法、半经验方法、DFT（Density Functional Theory,密度泛函数理论）方法，而经典力学模拟的方法主要依据分子力学、分子动力学、蒙地卡罗模拟、布朗动力学。量子力学是利用波函数来研究微观粒子的运动规律的一个物理学分支学科，它以分子中电子的
7、非定域化（delocalization）为基础，一切电子的行为以其波函数（wave function）表示。根据海森伯（Heisenberg）的测不准原理(uncertainly principle)，量子力学可计算区间内电子出现的概率，其概率正比于波函数绝对值的平方，而欲得到电子的波函数，则需解薛定谔方程式（Schrdinger equation）,即：（1. 1） = 此式中， 为薛定谔算子（是一些数学指令）、 为电子的波函数、E 为能 量。但是，分子中含有许多电子，要解此方程并非易事。今日，许多科学家们一方面致力于改进量子计算的方法及增进其精准度，如 1999 年诺贝尔物理学奖所颁给的密
8、度泛函数理论（density functional theory ,DFT）,即为非常精确的量子计算方法；另一方面，科学家们积极研究提升计算机的计算速度，以期计算较多电子的分子系统。通过量子力学，我们可以研究离子、原子、凝聚态物质，以及原子核和基本粒子的结构、性质，它与相对论一起构成了现代物理学的理论基础，也正是量子力学的出现，许多现象才得以真正被解释，新的、无法被直觉想象出来的现象被预言，而且基本上所有的物理间的基本相互作用都可以在量子力学的框架内描写。随着量子力学理论的逐步完善、经验力场的不断开发和更快更准确的算法以及计算机的普及和计算速度、计算量的不断提升，半个世纪以来，分子模拟的理论和
9、方法得到了快速的发展，在化学化工、材料科学、医药科学、生命科学等诸多领域发挥着越来越重要的作用，并形成了一门专门的学科分子模拟。今天，国内外的绝大多数的优秀科研机构和高等院校都已或正在涉足分子模拟领域，每年在世界范围内发表的相关研究论文数以万计，并且呈不断上升趋势。分子模拟所担当的角色也由早期纯粹的解释型逐渐过渡到解释、指导及预测并重型。通过上文，我们可以知道量子力学的方法适用于简单的分子或电子数目较少的体系。但是在自然界与工业上的许多系统中，比如聚合物（脂肪、橡胶、安全玻璃等等）、生化分子（核酸、酶、蛋白质、多糖类等）、合金材料等均含有大量的原子及电子。对于金属材料、聚合物材料、浓稠溶液、固
10、态混合物、纳米材料等系统，我们不需要了解单一分子的性质及分子间的相互作用，我们需要了解的是整个系统的各种集合的性质、动态行为与热力学的性质。对于这样的复杂体系，因其含有的电子数目较多，迄今为止，仍不可能完全仰赖量子力学计算。为解决这样复杂庞大的体系，科学家们在二十世纪六十年代就开始着手研究各种可行的非量子力学计算方法，利用这些非量子力学的计算方法可以解决许多复杂系统的问题，并得到相当精确的结果。而这些非量子力学计算方法通常指的就是经典力学模拟方法，即：分子力学、分子动力学、蒙地卡罗模拟、布朗动力学。分子力学方法（molecular mechanics ,MM）是依据经典力学（classical
11、 mechanics）的计算方法。因为用量子力学的方法去模拟计算大体系，如上千个原子的聚合物，是不实际的，即便是这种模拟能够进行，大多数情况下得到的大部分信息都要被舍弃。这是因为在模拟大体系的时候，常常是为了得到统计性质，如模量、扩散系数等，而这些量主要依赖于原子核的位置，此时计算量极大的电子运动信息并不重要。根据量子力学的波恩-奥本海默近似，分子的能量可以近似看做构成分子的各个原子的空间坐标的函数，简单的讲就是分子的能量随分子构型、相对位置变化而变化，而描述这种分子能量和分子结构之间关系的就是分子力场函数（force field）。而分子力学就是主要依据分子的力场计算分子的各种特性。力场函数
12、中含有许多参数，这些参数可由量子力学计算或通过科学实验获得。利用分子力学的方法可计算体系庞大的分子的稳定构象、热力学特性及振动光谱等信息。与量子力学相比较，此方法简便、快速，往往可快速得到分子的各种性质。在某些情况下，由分子力学方法所得到的结果几乎可以与高阶量子力学方法所得到的结果一致，但其所需的计算时间却远远小于量子力学的计算。故分子力学方法常被用于药物、团簇体、聚合物大分子的研究。分子动力学模拟（molecular dynamics simulation ,MDS）是目前最广泛采用的计算庞大复杂体系的方法。自 1966 年起，由于分子力学的发展，人们又系统的建立了许多适用于生化分子体系、聚
13、合物、金属及非金属材料的力场，使得计算复杂体系的结构与一些热力学、光谱性质的能力及精准性大为提升。分子动力学模拟是应用这些力场及根据牛顿运动力学原理所发展的计算方法。此方法的优点为精准性高，同时可以获得系统的动态与热力学统计资料，并可广泛地适用于各种系统及各类特性的研究。与蒙地卡罗计算方法相比较，分子动力学模拟系统中粒子的运动有正确的物理依据，并且计算技巧经过许多修改，现已日趋成熟，由于其计算能力强，能满足各类问题的需求，许多使用方便的分子动力学模拟商业软件也已问世。在先进国家的学校、工厂、医院等的实验室里，这些商业软件是不可或缺的研究工具。虽然分子动力学模拟具有很多优点，但是由于分子动力学模
14、拟计算需要引用数理积分方法，因此仅能研究系统在短时间内的运动，无法模拟一些运动时间长的运动问题（如由很多氨基酸组成的大分子蛋白质的折叠）。蒙地卡罗计算方法（Monte Carlo method）为最早针对庞大系统所采用的非量子计算方法,即 MC 计算法。蒙地卡罗计算法借由系统中质点（原子或分子）的随机运动，结合统计力学的概率分配原理来得到体系的统计及热力学资料，即根据待求问题的变化规律，构造合适的概率模型，然后进行大量统计实验，使模型的某些统计参量正好是待求问题的解。Monte Carlo 模型的建立可以分为三个具体的操作步骤，首先将所研究的物理问题演变为类似的概率或统计模型；其次通过数值随机
15、抽样实验对概率模型进行求解，其中包括大量的算术运算和逻辑操作；最后使用统计方法对得到的结果进行分析处理。这种计算方法可以研究复杂体系及金属的结构及其相变化性质。蒙地卡罗计算的弱点在于只能计算统计平均值，无法得到系统的动态信息。此计算所依据的随机运动并不适于物理学的运动原理，且与其他的非量子计算方法相比并非特别的经济快速，因此自从分子动力学计算逐渐盛行后，使用蒙地卡罗计算方法的人已逐渐减少。还有一种与分子动力学模拟类似的计算方法为布朗动力模拟（Brownian dynamics simulation ,BD）。布朗动力模拟适用于大分子的溶液系统，计算中，将大分子的运动分为根据力场作用的运动与来自
16、溶剂分子的随机力作用。通过布朗运动方程式可得到大分子运动的轨迹及一些统计与热力学的性质。布朗动力模拟通常适用于计算生化分子（如 DNA、RNA 、蛋白质、多肽）的水溶液。此算法的优点在于能计算大分子在较长时间范围内（纳秒级）的运动，其缺点为将溶剂分子的运动视为布朗运动粒子的假设未必是正确的。第 2 章 分子动力学模拟法分子模拟作为一种计算机模拟技术，主要可以进行解释工作和预测工作。前者为实验奠定理论基础，通过模拟解释实验现象、建立理论模型、探讨过程机理等，后者为实验过程提供可能性和可行性研究，进行方案辅助设计、材料性能预测、过程优化筛选等。分子模拟所涉及的研究领域，涵盖了物理、化学、化工、材料
17、、生化等几乎一切可以通过建立理论模型进行研究的体系，多数能够得到与实验结果近似的计算结果，所以分子模拟己经逐渐成为与实验技术并重的强有力的研究手段。20 世界 90 年代以来，由于计算机科学技术的飞速发展，分子模拟的地位日显突出。现在，理论分析、实验测定和模拟分析已成为现代科学的三种主要方法，在化工新产品、新材料和新药物的研究和开发中，采用分子模拟技术，可以从分子的微观性质推算及预测产品及材料的介观、宏观性质，已成为新兴学术方向。为探求纳米尺度下传热现象的规律和内在机理，我们必须从微观细节着手，研究载热粒子的行为，并依据统计力学原理得到系统的统计性质。分子模拟技术为某些问题( 如统计力学)的精
18、确解答起到非常重要的作用，否则这些问题很难或只能近似求解。此外，一些在极温极压条件下很难或者根本无法进行的实验，利用分子模拟却比较容易做到，在这种情况下，它实际上是对理论应用范围的一个较好的检测。通过对比分子模拟的结果与真实的实验结果，分子模拟的结果可以检验理论的正确性，相反，真实的实验结果与理论可以检测模型构建的合理性，此两者相互联系相互进步，通过对比，我们可以修正所建立的模型，并从模拟结果中进一步加深对所研究问题的理解，有助于解释出现的新现象。2.1 分子动力学概述分子动力学模拟（Molecular Dynamics Simulation ,MDS）是目前最广泛计算庞大复杂体系所采用的方法
19、，并且现已广泛应用于计算物理、化学、材料、生物科学等领域内的一种计算机模拟技术，它的思想是将组成系统的微观粒子视为牛顿经典粒子，将所研究的系统视为经典多体系统，通过选择恰当的场函数来描述系统内微观粒子间相互作用的势函数及系统外加约束条件后，利用牛顿运动定律来求解微观粒子的牛顿运动方程。在 MDS 中，势能模型的选取是非常重要的环节。模拟能否成功取决于能否正确地选择势能模型。对于简单分子，常采用硬球、软球、Lennard-Jones 模型、Kihara 模型等等。对复杂分子，可以采用多中心的位置位置相互作用模型，但中心间的相互作用仍为简单模型，对带电分子或离子，还要引入库仑作用，常用的多体作用势
20、有 Stillinger-Weber(SW)作用势、 Tersoff 势等。自 1966 年起，由于分子力学的发展，人们又系统的建立了许多适用于生化分子体系、聚合物、金属及非金属材料的力场，使得计算复杂体系的结构与一些热力学、光谱性质的能力及精准性大为提升。分子动力学模拟是应用这些力场及根据牛顿运动力学原理所发展的计算方法。此方法的优点为精准性高，同时可以获得系统的动态与热力学统计资料，并可广泛地适用于各种系统及各类特性的研究。与蒙地卡罗计算方法相比较，分子动力学模拟系统中粒子的运动有正确的物理依据，并且计算技巧经过许多修改，现已日趋成熟，由于其计算能力强，能满足各类问题的需求，许多使用方便的
21、分子动力学模拟商业软件也已问世。2.2 分子动力学模拟的基本步骤分子动力学模拟（Molecular Dynamics Simulation ,MDS）的出发点是将组成系统的微观粒子视为牛顿经典粒子，将所研究的系统视为经典多体系统，通过选择恰当的场函数来描述系统内微观粒子间相互作用的势函数及系统外加约束条件后，利用牛顿运动定律来求解微观粒子的牛顿运动方程，的到系统随时间演进的微观过程，最后基于统计力学计算得到系统的各种参数。分子动力学模拟一般由下面几个步骤组成：1.选取要研究的系统及其边界和系统内粒子间的作用势（力场函数）；2.设定系统中粒子的初始位置和动量初始化；3.建立模拟算法，计算粒子间作
22、用力及各粒子的速度和位置，并求解Newton 运动方程，得到下一时刻粒子的空间位置和速度；4.当体系达到平衡后，依据相关的统计公式，获得各宏观参数和输运性质；2.2.1 选择恰当系统与边界分子动力学模拟是在一定的系综下进行的，常见的系综有微正则系综、正则系综、等温等压系综和等温等焓系综等。在统计物理中，系综（Ensemble）代表一群类似的体系的集合。对一类相同性质的体系，其微观状态（比如每个粒子的位置和速度）仍然可以大不相同。实际上，对于一个宏观体系，所有可能的微观状态数是天文数字。统计物理的一个基本假设（各态历经假设）是：对于一个处于平衡的体系，物理量的时间平均等于对应系综里所有体系进行平
23、均的结果。体系的平衡态的物理性质可以对不同的微观状态求和来得到。系综的概念是由约西亚威拉德吉布斯（J. Willard Gibbs）在 1878 年提出的。在分子动力学模拟过程中，经常采用的系综有微正则系综(Microcanonical Ensemble)、正则系综(Canonical Ensemble)。正则系综 (Canonical Ensemble)也称为 NVT 系综，它是统计力学中系综的一种。它代表了许多具有相同温度的体系的集合，该系统具有恒定的粒子数N、体积 V 和温度 T。正则系综是最普遍应用的系综。通常，系综内每个体系的粒子数和体积都是相同的。但每个体系都可以和系综内其他体系交
24、换能量。同时系综里所有体系的能量总和，以及所有体系的总个数是固定的。但是正则系综的恒温条件是不能自发实现的，因此在分子动力学模拟过程中经常会采用一定的措施来保证系统温度的恒定。在分子动力学模拟过程中，系统的温度 T又 Boltzmann 能量均分定理得出 i，即:（2. 1）=123 式中，N 代表区域中原子的数量，v i 为粒子 i 的速度， kB=1.3810-23 JK-1。此公式成立的前提是系统区域内达到热力学平衡，若系统区域内无法达到热力学平衡，温度的定义就失去了意义。对于 NVT 系综，调节温度恒定的有效方法是直接调节系统内微观粒子的速度。在这些条件下，当系综内所有体系被分配到不同
25、的微观态上，我们发现：每个微观态上的体系个数正比于 exp(-E)。其中 =1/kBT ，kB 是玻尔兹曼常数，T 是绝对温度。正则系综的配分函数是：（2. 2）（）=(-E)配分函数的对数就是亥姆霍兹自由能（Helmholtz Free Energy，符号 F），当配分函数计算出以后，平均能量可以直接从 ln Z 对 一阶导数中求得。（2. 3）=- （2. 4）=-Z微正则系综(Microcanonical Ensemble)又称为 NVE 系综，该系综由许多具有相同能量 E，粒子数 N，体积 V 的体系的集合。微正则系综可以通过平衡态分子动力学由系统内微观粒子的牛顿运动方程自动实验，比较
26、简单方便，因此多数平衡态分子动力学模拟是在微正则系综下进行的。对于分子动力学模拟除了要选定合适的系综之外，还需要对模拟体系加以适当的边界条件。这是因为材料的宏观性质是由大量粒子的相互作用所表现出来的，在模拟过程中，如果所建立的模型含有足够多的粒子，那么才能真实的模拟出宏观材料的性质。但是计算机的计算能力是有限的，分子动力学模拟时粒子数目一般不超过 10000 个，加上系综的某些限制性条件（例如 NVT 系综中，调节系综温度恒定的某些方法），模拟的粒子数目可能要比这个数还要小。这样少的粒子所组成的系统与真实的宏观材料体系是具有一定误差的，同时粒子数目少，其表面原子与体内原子之比显然很大，因此会造
27、成明显的表明效应。因此，为了减小小体系的影响且使计算中系统的密度维持恒定，在分子动力学模拟过程中通常会采用有效的边界条件。常用于分子动力学模拟中的边界条件有三种：周期性边界条件、对称性边界条件和固壁边界条件 i。在周期性边界条件下，位于图中央的盒子表示所模拟的系统，其周围的盒子与模拟系统具有相同的排列及运动，我们称之为周期性镜像（Periodic mirror image）系统，当计算系统中任意微粒移出盒子后，则必有一粒子由相对方向移入，这样的限制条件使得系统中的粒子数目维持恒定，密度不变，符合实际要求。在周期性边界条件下，计算系统中分子间的作用力变得复杂，我们通常采取最近镜像法（Neares
28、t mirror image）来简化问题。如 图 2- 1 周期性边界示意图 。图 2- 1 周期性边界示意图2.2.2 分子间作用力的描述力场分子力场（Forcefields ）根据量子力学的波恩-奥本海默近似可以将力场看作是粒子势能面的经验表达式，即将一个分子的势能近似看作构成分子的各个原子的空间坐标的函数，简单地讲就是分子的势能随分子构型的变化而变化，而描述这种分子势能和分子结构之间关系的就是分子力场函数。分子动力学模拟过程中计算系统内粒子间的作用力主要依靠的就是力场给各个粒子间分配的作用力模型，分子动力学模拟的主要依据是经典力学，那么分子动力学模拟结果的准确性也主要取决于力场函数对粒子
29、间相互作用力描述的正确与否。由于分子动力学模拟研究的广泛性，分子间相互作用模型并非完全相同，针对不同的分子类型、分子体系，分别有不同的力场函数可以使用。分子力学与分子动力学发展至今已有 40 多年，其计算的系统由原来最初的单原子分子系统延伸至多原子分子、聚合物分子、生化分子系统。计算时所需的力场也随着系统的复杂度的增加而增加，在执行分子动力学计算时，选择合适的力场十分重要。常见的力场有CVFF、 CFF91、CFF95、PCFF、AMBER 、CHARMM、COMPASS 和 MM 形态力场等。MM 形态力场为 Allinger 等所发展，此力场适合用于模拟有机化合物、自由基、离子等随着分子模
30、拟的进步 MM 形态力场也发展出 MM2、MM3 和MM4 形态力场等，MM 力场将一些常见的原子细分（例如，将碳原子分为sp3、sp 2、sp、酮基碳、碳自由基和碳阳离子等并且这些不同形态的碳原子具有不同的力场参数），并为这些原子补充更多的力场参数，以更精确详细的来描述相同类型、不同电子结构原子的受力。CVFF 力场其全名为一致性价力场（Consistent valence force field）,为Dauber Osguthope 等所发展。此力场以计算系统的结构与结合能最为准确，亦可以提供合理的构型能和振动频率，常用于计算氨基酸、水及含各种官能团的分子体系，后经过不断改进，CVFF 力
31、场已适用于计算各种多肽、蛋白质与大量的有机分子。AMBER 力场（ Assisted model building with energy minimization）为美国加州大学的 Peter Kollman 等所发展，此力场通常可以得到分子的几何结构、构型能、振动频率与溶剂化自由能（Solvation free energy）等数据,所以常用于较小的蛋白质、核酸、多糖等生化分子的模拟计算。CHARMM 力场（Chemistry at Harvard macromolecular mechanics）为美国哈佛大学所发展，此力场可研究许多分子系统，例如小的有机分子、溶液、聚合物、生化分子等。
32、此力场参数除来自计算结果与实验值的比对外，还引用了大量的量子计算结果。应用此力场能够得到与实验值相近的结构、作用能、构型能、转动能、振动频率、自由能与许多和时间有关的物理量。随着力场应用研究的深入，随后出现了第二代力场。第二代力场的形式远比经典里场复杂，并需要大量的力场参数。其设计的目的是为了更加精确的计算分子的各种性质、结构、光谱、热力学特性、晶体特性等，其力场参数除了引用大量的实验数据之外，并参照了精确的量子计算结果，尤其适用于有机分子系统。CFF91 、CFF95 、 COMPASS 与 PCFF 都属于第二代力场。一般情况下，分子的总能量为动能与势能的和，而分子的势能可表示为简单的几何
33、坐标的函数。复杂分子的总势能一般可分为各类型势能的总和，即由成键项（U bonded）、交叉项 (Ucross)和非键项(U nobonded)三部分构成，即：总势能=成键项势能 +交叉项势能 +非键项势能习惯上常用的数学表达式为：U=U bonded+Ucross+Unobonded成键项（Bonded term）包括键伸缩项(Ub)、键角弯曲项 (U)、二面角扭曲项(U)和离平面振动项(U )。交叉能量项(Cross energy term)主要指成键作用之间耦合引起的能量变化，包括键伸缩-键伸缩(Ubond-bond)、键伸缩-键角弯曲-键伸缩(Ubond-angle-bond)、键角弯
34、曲- 键角弯曲(Uangel-angle)、二面角扭转项 -键伸缩(Utorsion-bond)、二面角扭转项- 键角弯曲- 键角弯曲 (Utorsion-angle-angle)、键角弯曲-二面角扭转项-键角弯曲(Uangle-torsion-angle)、离平面振动-二面角扭转项-离平面振动(Uoop-torsion-oop)，一般情况下，交叉能量项对整体势能影响较小,所以在一般的力场中并不含有这些项,只有在计算结果精度要求较高时才有交叉能量项。交叉项通常形式比较复杂,计算起来非常耗时,尤其是在计算原子数目较多的体系时。非键项（Nobonded term）包括范德瓦耳斯作用（U vdw）、
35、库伦静电作用（U el）和氢键作用（ UH）等。以下着重强调几个重要的势能项：(1) 非键结范德瓦耳斯势能(U vdw) 若 A、B 两原子属于同一分子但其间隔多于两个连接的化学键（如 ACCB），或者两原子属于不同的分子，虽无直接连接成键但因为距离较近，则此两原子间有非键结范德瓦耳斯力（Nobonded van der waals force）。计算非键结作用时，通常将某个原子视为位于其原子核坐标的一点，力场中最常见的非键结势能形式为 Lennard-Jones(L-J)势能，这种势能又称为 12-6 势能，其数学表达式为：（2. 5）()=4()12()6式中 与 为势能参数（因原子对而异
36、），r 为原子对间的距离，r-12 项为排斥项，r-6 项为吸引项，当 r 很大时势能趋近于零，这表明原子对距离很远时已无非键结作用，如 图 2- 2 L-J 势能曲线， 表 2. 1 部分常见分子的 LJ 势能参数。 图 2- 2 L-J 势能曲线表 2. 1 部分常见分子的 LJ 势能参数原子或分子 （/k b）/ K r/ He 10.41 2.602Ar 141.6 3.350O2 126.3 3.382CH4 161.3 3.721C3H8 268.5 4.992注: 1 =0.1nm 。 (2)键伸缩势能 (Ub) 在分子中相互键结的原子形成化学键，如 CH，OH 键等，化学键的键
37、长并非恒定不变，而是在其平衡值附近呈小幅度的振动，这种振动所产生的势能称为键伸缩势能或键伸缩势能项（Bond stretching term）。键伸缩势能项的一般形式为(将此形式视为简谐振动)：（2. 6）=12(0)2式中 kb 为键伸缩的弹力常数，弹力常数越大，振动频率越高，振动越快；ri 和 分别表示第 i 个键的键长与第 i 个键的平衡键长。0表 2. 2 MM2 力场一些键的伸缩势能参数键结平衡键长 r/ 弹力常数kb/kcal/(mol )2Csp3Csp3 1.523 317Csp3Csp2 1.497 317Csp2=Csp2 1.337 690Csp3H 1.113 662C
38、sp2=O 1.438 777Csp3Nsp3 1.438 367注：1kcal=4.187 kJ, 1 =0.1nm 。 有些力场（如 MM3 力场）为了增加计算的准确性，除了简谐振动项外增加了非简谐振动项(除了二次的简谐振动项外，还增加了三次与四次的非间歇振动项)以更准确的描述键伸缩势能，如：（2. 7）=12(0)2+3(0)3+4(0)4(3) 键角弯曲势能(U ) 若分子中连续键结的三原子形成一定的键角。与键的伸缩一样，这些键角并非维持不变，而是在其平衡值附近呈小幅度的振荡（一般视为简谐振荡），通过键角弯曲项（Angle bending term）可以描述这项势能项，键角弯曲项的一般
39、形式为：（2. 8）=12(0)2式中 i 及 分别表示第 i 个键的键角及其平衡键角，k 为键角弯曲的弹力0常数。为了增加计算的精确性，我们可以像键伸缩势能项一样增加高次非简谐项以提高计算的精确度。(4) 二面角扭曲势能(U) 若分子中连续键结四个原子并形成二面角（Dihedral angle），如丙烷 CH3CH2CH3 中的 HCCC、H CCH 等，一般分子中的二面角易于扭动，所以需要一个势能项来描述这种势能，我们称之为二面角扭曲项（Torsion angle term），二面角扭曲项的一般形式为：（2. 9）=121(1+cos)+2(1-cos2)+3(1+cos3)式中， 为二面
40、角角度，V 1、V 2 和 V3 均为二面角扭曲项的弹力常数。(5) 离平面振动势能(U) 在分子中，有些原子有共平面的倾向，如丙酮CH3C(O)CH3 中，碳原子与氧原子的平衡位置位于共同平面，其他如苯环上的碳氢原子、烯类中的碳氢原子等均有共平面的倾向，通常共平面的四个原子的中心原子会离开平面小幅度的振动，所以我们使用离平面振动项（Out-of-plane bending term）来描述这种势能，离平面振动项的一般形式为：(2. 10)=12 2式中， 为离开平面振动的角度，k 为离开平面振动项的弹力常数。(6) 库伦静电势能（U el） 离子或带有部分电荷的分子，则这些带电荷的粒子间因相
41、对距离存在着静电吸引或排斥作用，我们通过库伦作用项（Columbic interaction term）来描述这种静电作用，库伦作用项的一般形式为：(2. 11)=,式中，q i 和 qj 为第 i 个粒子与第 j 个粒子所带的电荷，r ij 为所带电荷粒子间的距离，D 为有效介电常数（effective dielectric constant）。2.2.3 Lennard-Jones(L-J)作用势与 Tersoff 作用势对于单壁碳纳米的分子动力学模拟，由于单壁碳纳米管模型中既有两体相互作用也有三体相互作用，所以针对此模拟体系，采用了两种作用势，分别是Lennard-Jones 作用势与
42、Tersoff 作用势。Lennard-Jones(L-J)作用势（又称 L-J 势能函数或 6-12 势能函数）是计算化学中用来模拟两分子间作用势能的一个函数，是一种两体作用势。最早由曼彻斯特大学的数学家 John Lennard-Jones 于 1931 年提出。由于其解析形式简单所以被广泛使用，L-J 作用势的表达式为： （2. 12）()=4()12()6L-J 作用势相应的两体作用力为：（2. 13）（）=-()=-()=4(121213667)式中 为 L-J 作用势的势能阱的参数， 是 L-J 作用势的互相作用势能正好为零时的两体距离（也称为平衡常数）,r ij 为粒子 i、j
43、间的距离。Tersoff 作用势是一种典型的三体作用势，它已被广泛用于描述碳族化合物的原子间的作用力 i。Tersoff 作用势的公式如下：（2. 14）=()()-()其中：（2. 15）()= -1- 2()（2. 16）()= -1- 2/()（2. 17）()=1,+ （2. 18）=+2 ,=1+( (,)() ()（2. 19）()=1+22 22+(-cos)2式中，、 为 L-J 势能参数， ijk 表示粒子 i-j 键和粒子 i-k 键之间的夹角，键序 bij 包含多体相互作用。公式中参数的选取如表：表 2. 3 Tersoff 势函数参数eV Snm-1a n c d hR
44、nmDnm5.16441.4471.576919.6401.572410-70.727511/(2n)3.80494.3484-0.570580.1950.015通过 Tersoff 作用势函数就可对势函数求导来计算各个原子的受力情况了。2.2.4 公式的无量纲化在分子动力学模拟计算系统中原子或分子的运动时，会涉及到大量的浮点运算和指数运算，如原子质量以 g 为单位时，则原子质量的量纲为 10-22g；位置若以 cm 为单位时，其量纲为 10-8cm；积分步程以秒（ s）为单位时，其量纲为 10-1310-16 s。这些量的量纲都非常小，在计算机计算时会增加大量的浮点运算，可能在模拟中导致计算
45、的误差，因此在实际执行分子动力学模拟计算时都会对计算公式进行无量钢化处理（也可看作是简化单位）以减少误差。本文中涉及到的简化单位转换式有：能量： （2. 20）=时间： （2. 21）=温度： （2. 22）=原子间距离： （2. 23）=速率： （2. 24）=作用力： （2. 25）=左上角用*标记的量就是简化单位后的参数，在分子动力学模拟计算中，各变量应按照简化单位后的数值进行运算。2.2.5 粒子的位置和动量初始化在分子动力学模拟过程开始前，必须先确定所有粒子的初始位置和速度，为了能使模拟系统尽快达到平衡并减少计算所需的时间，分子的初始速度应接近真实情况。由于系统经过一定的模拟时间达平
46、衡后，分子的速度的分布近似满足 Maxwell 统计速度分布，所以采取 Maxwell 统计分布来赋予分子的初始速度是非常理想的，Maxwell 统计分布如下：（2. 26）()=( 2)1/2-1222该式提供了一个质量为 mi 的原子在温度 T 时，在 X 方向的速度 vix 的高斯分布。随着模拟时间的推进，模拟的系统将遵循能量最小化的原则达到稳定状态。2.2.6 求解粒子的 Newton 运动方程在分子动力学模拟中，必须求解式的牛顿运动方程以计算速度与位置。而运动方程的迭代方法及离散方法很多，其中以 Verlet 算法、Gear 校正预测法算法（Predictor-Corrector method） i最为常用。在 Verlet 方法中，将粒子的位置以泰勒公式展开，即：(2. 27)(+)=()+()+12!22()()2+将式中得 换成 - 可得： (2. 28)(-t )=()()+12!22()()2+将式 2.27 与式 2.28 相加得：(2. 29)(+)=-r(-t )+2()+22()()2因为 ,所以根据式 2.29，可由 t 及 t- t 的位置预测 t+t 时的22()=()位置。将式 1.7 与式 2.27 相减的速度式：(2. 30)()=12(+)(-t)此式表示时间 t 时的速度可由 t+t 与 t-t 的位置
