第6章-分子动力学方法

..页脚第6章分子动力学方法经典分子动力学方法无疑是材料，尤其是大分子体系和大体系模拟有效的方法之一。分子动力学可以用于NPT，NVE，NVT等不同系综的计算，是一种基于牛顿力学确定论的热力学计算方法。与蒙特卡罗法相比在宏观性质计算上具有更高的准确度和有效性，可以广泛应用于物理，化学，生物，材料，医学等各个领域。本章在介绍分子动力学的基本概念的基础上，简单介绍了分子动力学的基本思想，势函数分类和基本方程。然后介绍了分子动力学的常用系综和典型的NPT，NVE，NVT系综基本方程。结合材料建模中的基本简化方法和技巧，阐述了边界条件和时间积分的数值处理技巧。最后，利用统计力学的基本概念给出分子动力学的计算信息的解析方式。并且结合MaterialsExplore软件计算分析了CNT的几何结构稳定性。6.1引言分子动力学方法(MolecularDynamics,MD)方法是一种按该体系部的禀动力学规律来计算并确定位形的变化的确定性模拟方法。首先需要在给定的外界条件下建立一组粒子的运动方程，然后通过直接对系统中的一个个粒子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计力学方法得到多体系统的静态和动态特性，从而获得系统的宏观性质。可以看出，分子动力学方法中不存在任何随机因素，这个也是分子动力学方法和后文要提到的蒙特卡洛方法的区别之一。在分子动力学方法的处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学定律(或者是拉格朗日方程)。每个分子运动的禀动力学是用理论力学上的哈密顿量或者拉格朗日函数来描述，也可以直接用牛顿运动方程来描述。确定性方法是实现玻尔兹曼的统计力学途径。这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。但是分子动力学方法的计算机程序相对蒙特卡罗较复杂，其计算成本较高。分子动力学方法发展历史改革经历了近60年，分子动力学方法是20世纪50年代后期由AlderBJ和WainwrightTE创造发展的。Alder和Wainwright在1957年利用分子动力学模拟，发现了“刚性球组成的集合系统会发生由其液相到结晶相的相转变”，后来人们称这种相变为Alder相变。其结果表明，不具有引力的粒子系统也具有凝聚态。到20世纪70年代，产生了刚性体系的动力学方法，成功地被应用于水和氮等分子性溶液体系的处理；1972年，LessAW和EdwardsSF等发展了该方法并扩展到了存在速度梯度，即处于非平衡状态的系统。之后，此方法被GinanMJ等推广到了具有温度梯度的非平衡系统，从而构造并形成了所谓的非平衡分子动力学方法体系。进人20世纪80年代之后，出现了在分子部对一部分自由度施加约束条件的新的分子动力学方法，从而使分子动力学方法可适用于类似蛋白质等生物大分子的解析与设计。分子动力学方法真正作为材料科学领域的一个重要研究方法，开始于由Andersen,Parrinello和Rahman等创立恒压分子动力学方法和Nǒse等完成恒温分子动力学方法的建立及在应用方面的成功。后来，针对势函数模型化比较困难的半导体和金属等，1985年人们又提出了将基于密度泛函理论的电子论和分子动力学方法有机统一起来的所谓Car-Parrinello方法，亦即第一性原理分子动力学方法。这样，分子动力学的方法进一步得到发展和完善，它不仅可以处理半导体和金属的间题，同时还可应用..页脚于处理有机物和化学反应。关于Car-Parrinello分子动力学方法将在第7章重点学习讨论。本章重点讨论经典分子动力学方法的基本原理和计算方法。6.2分子动力学的计算框架6.2.1基本思想分子动力学方法是沿用牛顿运动方程或拉格朗日方程、哈密顿方程，通过考察粒子的运动来研究多粒子系统的物理性质。在处理孤立粒子或原子团簇(不考虑与其他粒子的相互作用)时，单纯地对牛顿运动方程进行积分求解就行了；而在处理凝聚系统时，可以认为将所考虑的N个粒子放入具有一定体积的容器中，在这样组成的封闭系中，建立直角坐标(x,y,z)，每个粒子的位置就由三个坐标分量决定，通过求解3N个联立方程组就可以得到保守系的总能量。分子动力学的基本思想在于通过设定原子之间的相互作用(势函数)和相关的系综(亦即作用对象和条件)，确定其基本的模拟畴。在给定的牛顿方程、拉格朗日方程或者是哈密顿方程进行时间的迭代，在达到指定的力学收敛条件之后得到一个最终的位置坐标，然后通过相关的位形信息和运动的速度和加速度等信息通过热力统计学的方法，给出模拟体系的统计信息，主要包括复杂的光学，电学和一些其他的物理信息。其基本思想，总结如图6-1所示。势函数，相互作用系综，外界条件运动方程（牛顿方程，哈密顿方程，拉格朗日方程）原子的坐标位置原子的坐标和速度热力学性质动力学性质光学性质原子运动三维结构一次信息二次信息动力学方程初始化条件图6-1分子动力学的一般思想6.2.2计算流程总体来说，分子动力学的基本计算有如下四个步骤：1.模型的设定，也就是势函数的选取势函数的研究和物理系统中对物质的描述研究息息相关。分子动力学的模拟最早使用是硬球势，即小于临界值时无穷大，大于等于临界值时为零。常用的是Lennard-Jones势函数，还有嵌入原子EAM势函数，不同的物质状态描述用不同的势函数。模型势函数一旦确定，就可以根据物理学规律求得模拟中的守恒量。关于势函数选取将在6.4节进一步讨论。2.给定初始条件运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。如Verlet算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步的坐标或者一组零时刻的速度值。..页脚一般意义上讲系统的初始条件不可能知道，实际上也不需要精确选择代求系统的初始条件，因为当模拟时间足够长时，系统就会忘掉初始条件(对于无记忆的体系而言)。当然，合理的初始条件可以加快系统趋于平衡的时间和步程，获得好的精度。常用的初始条件可以选择为，令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到；令初始位置随机的偏离差分划分网格的格子上，初始速度为零；令初始位置随机的偏离差分划分网格的格子上，初始速度也是从玻尔兹曼分布随机抽样得到。启动计算设定坐标和相互作用t+Δt计算作用在原子上的力计算物理量并对其结果进行统计结束ttmaxYN(/2)(/2)()/()()(/2)iiiiiiiVttVtttFtmRttRttVtt图6-2分子动力学计算流程3.趋于平衡计算在边界条件和初始条件给定后就可以解运动方程，进行分子动力学模拟。但这样计算出的系统是不会具有所要求的系统能量，并且这个状态本身也还不是一个平衡态。为使得系统平衡，模拟中设计一个趋衡过程，即在这个过程中，增加或者从系统中移出能量，直到持续给出确定的能量值，称此时系统已经达到平衡，达到平衡的时间称为弛豫时间。分子动力学中，积分步程的大小选择十分重要，这决定了模拟所需要的时间。为了减小误差，步长要小，但过小系统模拟的弛豫时间就长。因此根据经验选择适当的步长。如，对一个具有几百个氩气分子的体系，采用Lennard-Jones势函数，发现取h为0.01量级，可以得到很好的结果。这里选择的h是没有量纲的，实际上这样选择的h对应的时间在10-14s的量级。如果模拟1000步，系统达到平衡，弛豫时间只有10-11s。4.宏观物理量的计算实际计算宏观的物理量往往是在模拟的最后阶段进行的，它是沿相空间轨迹求平均来计算得到的，根据各态历经假说，时间平均代替系综平均。关于各态历经假说将在第8章的蒙特卡洛方法中深入讨论。..页脚6.2.3经典分子动力学中近似处理实际的计算体系由于其外在条件的复杂性，往往需要对计算的系统进行一定的近似处理，在经典分子动力学计算中，引进以下近似处理：经典粒子相互作用，不考虑电子相互作用量子效应。力的作用形式，由参数可调的相互作用势函数决定，并经过实验测定来验证。模拟体系与实际体系相差较大，一般需要采用周期边界来扩展计算体系。时间平均是在有限时间完成。【练习与思考】6-1.查找文献，根据上述的分子动力学的处理流程图，编写实现分子动力学的简单程序，可参考DaanF和BerendS编著的《分子模拟》一书。6.3分子动力学的系综系综(ensemble)是指在一定的宏观条件下(约束条件)，大量性质和结构完全相同的、处于各种运动状态的、各自独立的系统集合，或称为统计系综。系综是用统计方法描述热力学系统的统计规律性时引入的一个基本概念；系综是统计理论的一种表述方式，系综理论使统计物理成为普遍的微观统计理论；系综并不是实际的物体，构成系综的系统才是实际物体。约束条件是由一组外加宏观参量来表示。在平衡统计力学畴下，可以用来处理稳定系综。6.3.1常用系综分类1.正则系综(canonicalensemble)全称应为“宏观正则系综”，简写为NVT，即表示具有确定的粒子数(N)、体积(V)、温度(T)。正则系综是分子动力学方法模拟处理的典型代表。假定N个粒子处在体积为V的盒子，将其埋入温度恒为T的热浴中。此时，总能量(E)和系统压强(P)可能在某一平均值附近起伏变化。平衡体系代表封闭系统，与大热源热接触平衡的恒温系统。正则系综的特征函数是亥姆霍兹自由能(,,)FNVT。2.微正则系综(micro-canonicalensemble)简写为NVE，即表示具有确定的粒子数(N)、体积(V)、总能量(E)。微正则系综广泛被应用在分子动力学模拟中。假定N个粒子处在体积为V的盒子，并固定总能量(E)。此时，系综的温度(T)和系统压强(P)可能在某一平均值附近起伏变化。平衡体系为孤立系统，与外界即无能量交换，也无粒子交换。微正则系综的特征函数是熵(,,)SNVE。3.等温等压(constant-pressure,constant-temperature)简写为NPT，即表示具有确定的粒子数(N)、压强(P)、温度(T)。其总能量(E)和系统体积(V)可能存在起伏。体系是可移动系统壁情况下的恒温热浴。特征函数是吉布斯自由能(,,)GNPT。4.等压等焓(constant-pressure,constant-enthalpy)..页脚简写为NPH，即表示具有确定的粒子数(N)、压强(P)、焓(H)。由于HEPV，故在该系综下进行模拟时要保持压力与焓值为固定，其调节技术的实现也有一定的难度，这种系综在实际的分子动力学模拟中已经很少遇到了。5.巨正则系综(grandcanonicalensemble)简写为μVT，即表示具有确定的粒体积(V)、温度(T)和化学势(μ)。巨正则系综通常是蒙特卡罗模拟的对象和手段。此时、系统能量(E)、压强(P)和粒子数(N)会在某一平均值附近有一个起伏。体系是一个开放系统，与大热源大粒子源热接触平衡而具有恒定的温度(T)。特征函数是马休(Massieu)函数J(μ，V，T)。对于不同的系综，由于其选择的运动方程的求解方法亦有不同。针对不同的计算体系具体的细节比较复杂，下文将以NEV和NTP这两类简单的系综运用两类不同的条件方法下求解的不同方程做简要的介绍。6.3.2NEV系综基本方程NEV系综是一类较为常见的系综。对于该系综，经典分子动力学方法一般采用差分近似法求解牛顿运动方程，并追踪系统的时间变化。考虑由N个粒子构成的系统，假设第i个粒子在时刻()ixt位置，速度为()ivt受到的作用力为()iFt，为求出该粒子(i)在时刻tt的位置，经常使用的方法是Verlet算法(对于该算法的具体细节在下文将有详细的推导)。根据系统的动能及统计物理学中的能量均分定理，则系统的平衡温度可表达为2212iiBfipTkNm(6.1)根据此方法，在晶体中粒子的配置与排布、晶形变化等结构相变也可以用计算机模拟进行研究。根据物理化学的基本气体的基本方程可以知道，体系的NEV系综下的计算体系主要变化的量有两个压力和温度。下文就从恒