分子模拟作为理论与计算化学的一个重要方法，可以从分子微观层面上认识和解释各种物理化学现象、化学反应机理和分子的物理化学性质，结合统计力学对宏观体系的热力学和动力学性质进行模拟和预测，逐渐成为人们认识微观世界的重要手段。量子力学方法尽管为我们提供了电子态的信息和准确方法，但由于工作量巨大限制了其应用的分子和体系尺度，对于生物大分子、溶液体系、超分子材料等复杂体系显然无能为力。经典分子力场方法忽略了电子的自由度，体系的能量只与原子核的位置有关，并假定体系中原子和分子的运动服从经典力学。分子势能简单分解为成键势能项和非键势能项，分子力场参数易从量子力学方法或实验方法得到，大大提高了模拟效率，分子力学方法为我们模拟生物大分子、凝固态和复杂体系的热力学和动力学性质提供了有力工具，在药物设计、生物分子结构与识别、材料设计、凝聚态性质等领域已经取得了富有成效的研究成果。经典分子力场模拟的可靠性关键取决于力场的质量，目前大部分常用力场，如OPLS，CHARMM，AMBER和GROMOS等，通常采用固定的原子中心固定点电荷来计算静电相互作用势能，无法体现这种极化影响，静电极化作用在气液界面、电解质溶液、生物分子溶液等凝聚态中发挥中至关重要的作用，发展极化力场已成为近年来力场发展的趋势之一。
本书详细介绍了分子力场的基本理论和基本方法，结合我们实验室多年来在极化分子力场方法方面的研究结果，重点阐述分子极化力场的分类、基本物理模型和相关的应用。简单阐述分子动力学模拟的基本原理和模拟细节，分类介绍水分子力场、离子水分子相互作用分子力场和蛋白质分子力场。结合和比较不同力场模型下的分子团簇、水溶液、电解质溶液和大分子构象等性质来具体阐述静电极化和力场模型的影响。在本书的写作过程中，张强主要负责第1章、第2章、第3章、第4章和第8章的写作，约24万字，张霞主要负责第5章、第6章和第7章的写作，约4.5万字。
特别感谢我们实验室的杨忠志和王长生老师，以及赵东霞、宫利东、吴阳、李欣等的无私帮助和付出的大量工作。同时感谢金振兴老师、赵岷老师、蔡克迪老师、张庆国老师在书籍出版过程中提供的建议和帮助。
分子模拟研究内容丰富，涉及面广，由于作者能力有限和时间仓促，书中难免有错误和不妥之处，敬请读者批评指正。
张强渤海大学2012年3月第1章分子模拟力场方法
1.1分子力学概述
1998年，美国化学家库恩（WalterKohn）和英国科学家波普尔（JohnA.Pople），因其对量子化学的贡献而获得诺贝尔化学奖，这对理论与计算化学科学来说具有里程碑意义，既是对其本人，也是对理论与计算化学领域的肯定，同时更是对该领域工作者和科学家的激励和鞭策。理论与计算化学作为化学的一个新兴分支，为整个化学的发展提供了强有力的理论支撑和科学指导。我们可以利用理论和计算化学来解释实验现象和结果，从分子微观层面上更加深入地了解各种化学现象、研究化学反应，对分子的物理和化学性质，乃至一些宏观体系的热力学和动力学性质进行模拟和预测，并通过数字可视化手段，让我们能更加直观地了解微观世界。
计算化学的首要方法就是量子力学方法[1，2]，量子力学是以波函数的形式描写分子中电子的非定域化(delocalization)行为，一切电子的行为以其波函数（wavefunction）表示，根据海森堡（Heisenberg）的测不准原理（uncertaintyprinciple），区间内电子出现的几率正比于波函数绝对值的平方，通过求解Schrdinger方程得到电子的波函数。最为普遍的量子力学计算方法为从头计算法(abinitio)。常用的分子轨道（molecularorbital）方法利用变分原理(variationprinciple)，将系统单电子的波函数展开为原子轨道波函数的组合，而原子轨道的波函数又是一些特定数学函数(如高斯函数)的组合。虽然这种方法计算精确，但工作量巨大，大约与电子数的五次方成正比，因而所能计算的体系受到了限制，通常不超过100个原子。为了提高量子力学计算的效率，自1960年起，陆续发展出一些近似的量子化学计算方法，这些方法多利用一些实验数据作为参数来取代真正的积分项，称之为半经验分子轨道计算方法（semi-empiricalmolecularorbitalmethod）。目前常用的半经验分子轨道计算方法有CNDO、INDO、MINDO、MNDO、ZINDO和AM1等[3，4]。利用半经验分子轨道方法可计算相对较大的分子，而且计算的结果往往可以与从头计算方法相媲美。但是，即使利用半经验方法和最先进的计算机技术（并行处理），目前所能计算的量子体系实际上也不超过1000个电子。科学家们正致力于改进量子计算方法及增进其精确度，密度泛函理论在近几年得到了快速的发展，在密度泛函理论中，对于电子体系的描述，使用电子密度(r)作为基本变量，而不像通常理论使用电子数目（N）和外势场（V）作为基本变量。电子密度决定着体系基态的波函数和体系的所有其他性质[5，6]。
简而言之，单电子密度(r)是决定体系一切性质的唯一自变量函数。以密度泛函理论为基础发展起来的计算体系能量和其他性质的方法，由于用三维的单电子密度代替了3N维的电子波函数，用求解比较简单的泛函方程代替了求解复杂的Schrdinger方程，计算工作量只约与电子数目（N）的三次方成比例，比通常量子化学从头计算（abinitio）方法所需要的时间大为减少。
1929年，Dirac指出：“大部分的物理和整个化学领域的根本定律以及相应的数学处理方法我们是清楚的，但是困难就在于这些方程太复杂了，很难求解。”量子力学的方法仅适用于简单的分子，或电子数目较少的体系。在下面几种情况下，利用量子力学方法来处理是很困难的。第一，生物大分子（蛋白质、酶、酵素、核酸、多糖类……）、聚合物（橡胶、安全玻璃、脂肪、油类分子……）等含有大量数目的原子及电子体系。第二，分子环境的影响。量子力学往往计算一些气态的分子，而大部分化学过程都发生在溶液或其他凝聚态体系中，对溶剂和环境的计算无疑是对量子力学方法的巨大挑战。第三，利用从头计算方法计算准确的热力学和动力学数据是很困难的，也只能局限于较小的体系，然而，这些结果可以帮助我们更好地指导工业生产和改进工艺。此外，金属材料、聚合物材料、浓稠溶液、纳米材料等系统，迄今仍不可能依赖量子力学计算。针对这样庞大的体系，科学家们又把目光投向了经典力学，从1960年左右就开始着手研究分子力学方法，以此来研究各种复杂体系，并取得了可喜的成果。
分子力学方法[7]（MolecularMechanics，MM）起源于1970年左右，是依据经典力学（ClassicalMechanics）的分子力场（ForceField）计算方法。分子力学的建立基于两条根本近似：第一，Born-Oppenheimer近似。由于组成分子体系的原子核的质量比电子大103~105倍，因而在分子中的电子运动速度将比原子核快得多，对于原子核位置的任何变化，迅速运动的电子都能立即进行调整，建立起与变化后的核力场相应的运动状态。从而，一方面可以把对原子和分子运动的数学描述同对电子结构及其运动的描述分离开来；另一方面，在分子内，源于分子电子结构变化的物理相互作用（如核与电子，电子同电子相互作用）都可以用依靠于组成分子的原子坐标的函数来描述。第二，体系中原子和分子的运动服从经典力学。也就是说，粒子的任何能级光谱都是连续的，而不是分立的；假定原子和分子的运动足够快，则可以忽略任何量子影响。认为体系的运动服从牛顿运动定律，而不是Schrdinger方程。
依据Born-Oppenheimer近似，计算中忽略电子的运动，将系统的能量视为原子核位置的函数。分子力场中有许多参数，这些参数可由量子力学方法或实验方法得到。
利用分子力学方法可以计算庞大和复杂分子的稳定构象、热力学性质、动力学性质及振动光谱等。与量子力学相比较，分子力学方法要简单得多，而且往往能够快速得到分子的各种性质。某些情形下，由分子力学方法所得到的结果几乎与高水平量子力学方法所得的结果一致，但其所需要的计算时间却远远小于量子力学的计算时间。目前，分子力学方法常被应用于药物、团簇体系和生物大分子等体系的研究。
最早对庞大体系采用非量子力学计算方法的是蒙特卡罗[8，9]（MonteCarloMethod），简称MC计算法。1953年，Metropolis和他的合作者报道了世界上第一例蒙特卡罗模拟，蒙特卡罗计算法是由系统中质点（原子和分子）的随机运动，结合统计力学的几率分配原理，得到体系的统计及热力学资料。MC计算方法至今仍经常被应用于研究复杂体系的结构及其相变化等性质。其弱点在于只能得到统计的平均值，无法计算系统的动态信息，而且它所依据的随机运动并不符合物理学的运动原理，与其他的非量子计算方法比较并非经济快速。但是，该方法能够相对高效地寻找体系低能构型。
分子动力学模拟（MolecularDynamicSimulation）方法[9，10]，简称MD计算法，是目前应用更加广泛的计算庞大复杂体系的方法。1957年，Alder和Wainwright报道了世界上第一例分子动力学模拟[11]；1964年，Rahman采用Lennard-Jones函数对于液体氩进行了分子动力学模拟；1971年，Stillinger和Rahman模拟了液态水的分子动力学过程[12]；1977年，McCammon等人报道了第一类蛋白质分子动力学模拟；1983年，Levitt完成了第一类核酸的分子动力学模拟[13，14]。自1970年起，由于分子力学的迅速发展，又系统地建立了许多适用于生化分子体系、聚合物、金属和非金属材料的力场，使得计算复杂体系的结构、热力学和光谱性质的能力及准确性大为提升。分子动力学模拟是应用这些力场及牛顿运动力学原理发展起来的计算方法，其优点在于系统中粒子的运动有正确的物理依据，精确性高，可同时获得系统的动力学和热力学统计资料，并可广泛地适用于各种系统及各类特性的探讨。
另有一种与分子动力学模拟类似的计算方法是布朗动力学模拟(BrownianDynamicsSimulation，BD)。布朗动力学模拟适用于大分子的溶液体系，计算中将大分子的运动分为依力场作用的运动和来自溶剂分子的随机力作用的运动。利用求解布朗运动方程式可得到大分子体系的运动轨迹及一些统计和热力学的性质。布朗动力学模拟通常用于计算生化分子(多肽、蛋白质和DNA等)的水溶液。此方法的优点在于能够在较长时间范围内（纳秒，ns）计算大分子体系的运动。其缺点是：将溶剂分子的运动视为布朗运动粒子，该假设缺乏合理性。
1.2分子力场作用项的一般形式
基于分子力学的基本近似，只要某一势能函数能够很好地描述体系内的各种相互作用，体系内的原子和分子的运动也可以用源于该势能函数的经典力来描述。那么，该势能函数以及相关的参数就称为力场。力场方法发展经历了近七八十年的历程，1930年，Andrews提出分子力学的基本思想[15]；1961年，Hendrickson完成了对大于六元环分子的构象分析[16]；1965年，Wiberg首次出版了能够搜寻分子能量最低点的通用力场程序；1976年，Allinger发表了MM1力场[17，18]；1981年，Kollman发表了第一个AMBER力场[25]；1983年，Karplus发布了分子力学与分子动力学模拟程序CHARMM[26]。分子力场不断完善，根据要研究的体系和性质的差异，更加具有针对性。经典分子力学模拟结果的优劣取决于所采用的分子力场的质量，所以发展更好的力场是分子力学经典模拟的关键。
事实上，每个分子都有其固定的力场，在实际应用时，将力场分解成不同的组分，用理论计算和实验拟合的方法建立力场参数，并且要求这套力场参数具有较强的通用性，可以对一类分子进行计算。力场是在光谱分析时为分析和预测振动光谱发展起来的，早期的分子力场是对振动力场的改进，侧重于计算分子的结构和能量，但不能预测分子的振动光谱。主要理论基础为：一个由n个原子组成的分子，坐标为X(x1，x2，，xn)，当分子的平衡构型变为其他构型时，势函数对其平衡坐标的泰勒展开为：
展开式中的第一项V0是分子在平衡时的势能，可设置为0；第二项在V取极小值时也为0；当分子构型变化不大时，二次项以上的各项也可忽略不计，这时势函数主要决定于第三项：
ijf称为力场参数。用ijf（i，j=1，2，3，，3n）可组成一个矩阵，如果忽略非对角项，则满足胡克定律，称分子体系为非耦合体系。如果分子的构型变化较大，必须考虑上式中的高阶项，此时分子体系处于非谐振态。1968年，Lifson和Warshel发展了自治力场(ConsistentForceField，CFF)[19]，目的是为了既能模拟分子的结构，又能模拟分子的振动光谱。他们用最小二乘法拟合理论计算和实验结果来获取力场参数，其为力场的参数化奠定了基础。
