5.分子动力学方法

5.分子动力学原理与方法2013年9月3日计算材料学本课提纲本节课将向同学位介绍分子动力学的原理和方法。分子动力学简介分子动力学的原本原理分子动力学的计算方法分子运动方程的数值求解边界条件与初值物质的势函数1.分子动力学简介分子动力学（MolecularDynamics，简称MD）是用计算机来模拟统计力学的一种计算方法。MD模拟用来研究不能用解析方法来解决的复合体系的平衡性质和力学性质，在数学、生物、化学、物理学、材料科学和计算机科学交叉学科占据重要地位。对MD使用者来说，牢固掌握经典统计力学、热力学系综、时间关联函数以及基本模拟技术是非常必要的。分子动力学模拟方法是一种确定性方法，是按照该体系内部的内禀动力学规律来确定位形的转变，跟踪系统中每个粒子的个体运动；然后根据统计物理规律，给出微观量(分子的坐标、速度）与宏观可观测量（温度、压力、比热容、弹性模量等）的关系来研究材料性能的一种方法。1.分子动力学简介分子动力学方法首先需要建立系统内一组分子的运动方程，通过求解所有粒子的运动方程，来研究该体系与微观量相关的基本过程。根据波恩-奥本海默近似，将电子的运动与原子核的运动分开处理，电子的运动用量子力学的方法处理，而原子核的运动用经典动力学方法处理。原子核的运动满足经典力学规律，用牛顿定律来描述，这对于大多数材料来说是一个很好的近似。对于一些较轻的原子和分子的平动、转动或振动频率n满hnkBT时，才需要考虑量子效应。1.分子动力学简介使用分子动力学方法计算时，需要三个步骤：建立一个由N个粒子（分子）组成的模型体系，包括系统中的原子，原子间的相互作用势，初始条件和边界条件；求解N个粒子(分子)组成的模型体系的牛顿运动方程，并用统计力学的方法记录每一步原子位置和速度;对模拟结果进行统计分析，得到材料的结构变化和各项性能，如力学性质，物理性能等。1.分子动力学简介Alder和Wainwright于1957年和1959年提出并应用于理想“硬球”液体模型，发现了由Kirkwood在1939年预言的“刚性球组成的集合系统会发生由液相到结晶相的转变”。Rahman于1963年采用连续势模型研究了液体的分子动力学模拟。Verlet于1967年给出了著名的Verlet算法，即在分子动力学模拟中对粒子运动的位移、速度和加速度的逐步计算法。1980年Anderson研究了恒压状态下的分子动力学，提出了等压分子动力学模型。1981年，Parrinello和Rahman给出了恒定压强的分子动力学模型，将等压分子动力学推广到元胞的形状可以随其中粒子运动而改变的范围。Nose于1984年提出了恒温分子动力学方法。1985年Car和Parrinello提出了将电子运动与原子核运动一起考虑的第一性原理分子动力学方法。2.分子动力学的原本原理分子动力学模拟是一种用来计算经典多体体系的平衡和传递性质的一种确定性方法。分子动力学中处理的体系的粒子的运动遵从牛顿方程：F=ma。原子i所受的力Fi可以直接用势能函数对坐标r的一阶导数，因此对N个粒子体系的每个粒子有这些方程的求解一般来说需要通过数值方法进行，这些数值解产生一系列的位置与速度对。要求解此方程组，必须要给出体系中的每个粒子的初始坐标和速度。2.分子动力学的原本原理给出原子的初始坐标和初始速度后，以后每一时刻的坐标和速度都可以确定。分子动力学整个运行过程中的坐标和速度称为轨迹（trajectory)。数值解普通微分方程的标准方法为有限差分法。给定t时刻的坐标和速度以及其他动力学信息，那么就可计算出t+Dt时刻的坐标和速度。Dt为时间步长，它与积分方法以及体系有关。输人信息包括原子间或分子间相互作用力的势函数，以及温度和压力等物理环境条件。2.分子动力学的原本原理对得到各时刻的原子位置坐标进行统计计算，则可得到有关的热力学性质(热力学能、比热容等）；同时统计处理各时刻的原子位置坐标和速度，就可得到动力学性质(扩散系数、粘滞系数等)。2.分子动力学的原本原理分子动力学计算的过程如下：①输入指定运算条件的参数(初始温度，粒子数，密度，时间步长）；②体系初始化(选定初始坐标和初始速度）；③计算作用在所有粒子上的力；④解牛顿运动方程(第③和第④步构成了模拟的核心，重复这两步，直到体系的演化到指定的时间。⑤经过统计力学，计算并输出物理量的平均值，及随时间的变化。2.分子动力学的原本原理分子动力学方法只考虑多体系统中原子核的运动，忽略电子运动的量子效应，在很宽的材料体系都较精确；对于涉及电荷重新分布的化学反应、键的形成与断裂、解离、极化以及金属离子的化学键都不适用，此时需要使用量子力学方法。经典分子动力学方法(MD)也不适用于低温，因为量子物理给出的离散能级之间的能隙比体系的热能大，体系被限制在一个或几个低能态中。当温虔升高或与运动相关的频率降低时，离散的能级描述变得不重要，更高的能级可以用热激发描述。对牛顿物理特征运动频率的初略估计可以根据简谐分析进行。对简谐振动来讲，量子化的能量为hn。显然，经典方法对相对的高频率的运动不适用。3.分子动力学的计算方法分子动力学计算的流程建立计算模型，设定计算模型的初始坐标和初始速度;选取合适的原子间相互作用势函数，以便进行力的计算；选择合适的算法、边界条件和外界条件;选定合适的时间步长，计算坐标和速度的改变；对计算数据进行统计处理。执行预定步数的循环计算;3.分子动力学的计算方法分子动力学由三个主要部分组成:初始化，运动轨迹和结果分析。初始化要求给每个粒子给定初始坐标和速度。即使初始坐标和速度可以从实验中得到，指定的开始矢量也不一定对应于所使用的势函数的最小值，因此需要进一步的最小化来弛豫应力。初始速度矢量根据伪随机数进行设置，使体系的总动能与目标温度对应，根据经典的能均分定理，在热平衡时，每个自由度的能量为kBT/2。可以通过给速度分量设置麦克斯韦分布或正态分布。3.分子动力学的计算方法当体系的初始坐标和初始速度设置以后，在进行模拟体系的性质以前，首先必须使体系进行趋于平衡的过程。在这个过程中，动能、势能相互转化，当动能、势能、总能量只在平均值附近涨落时，体系就达到平衡了。MD轨迹的无序性是经常出现的，无序行为意味着初始条件的一个微小的变化就会导致在很短时间内指数偏离正常轨迹。初始条件变化越大或时间步长越长，不稳定性发生的越快。在短时间内，偏差为指数增加，接下来是线性绕阈值振动。大多数MD模拟需要很长的时间步为几万到几十万时间步。3.分子动力学的计算方法时间步长和势函数时间步长Dt的选取是非常重要的，不合适的时间步长可能会导致模拟的失败或结果的错误，或者造成模拟的效率太低。时间步长的选取参考原子或分子特征运动频率。势函数是描述原子(分子）间相互作用的函数。原子间相互作用控制着原子间的相互作用行为，从根本上决定材料的所有性质，这种作用具体由势函数来描述。在分子动力学模拟中，势函数的选取对模拟的结果起着决定性的作用。在选取势函数时，不可盲目选取，要认真考察势函数是否合适。通过文献调研，查明它的出处、应用范围；还要自己验证势函数的好坏，通过模拟材料一些已知的性质来验证，然后再进行使用。3.分子动力学的计算方法力的计算方法在分子动力学模拟中所花费的计算时间，其中90%以上是用来计算作用在原子上的力，所用时间大致正比于原子数目的平方。一般分子动力学模拟的原子数目较大(几万、几十万甚至更大），因此对原子作用力进行简化求解是非常必要而有意义的。对于短程力，采用截断半径法，即只计算力程以内的作用力就可以了，这样大大减少了计算量；而对于像库仑力这样的长程力，需要找到近似处理办法来减少计算量，其中Ewald求和法就是常用的一种。3.分子动力学的计算方法截断半径法预先选定一个截断半径rc，只计算以截断半径为球体内的粒子间的作用力；与粒子之间的距离超过截断半径rc时，则不考虑它们的作用。选择截断半径rc时，L2rc(L为分子动力学的周期性盒子的长度）；当系统粒子数量很大时，还可以利用邻域列表法判断粒子的分布情况，进一步节省机时。对系统中的每个粒子建立邻域表。计算过程中，每隔一定步数需要更新此表。3.分子动力学的计算方法Ewald加和法Ewald加和法是Ewald于1921年在研究离子晶体的能量时发展的。一个粒子与所有其他原子发生作用：这个式子的求和收敛得非常慢，Ewald在求和时采用了一个技巧，加快了级数收敛的速度。3.分子动力学的计算方法算法的选取在分子动力学的发展过程中，人们发展了很多算法，常见的方法有Verlet算法、Leap-frog和VelocityVerlet算法、Gear算法、Tucterman和Berne多时间步长算法。在这些算法中，选取的准则、其评判标准中最重要的一点是能量守恒。动能和势能的涨落总是大小相等，符号相反。RMS涨落大约为0.025J/mol，动能和势能的RMS涨落大约为10.45J/mol。4.分子运动方程的数值求解多粒子体系的牛顿方程无法求解析解，需要通过数值积分方法求解，运动方程可以采用有限差分法来求解。将积分分成很多小步，每一小步的时间固定为dt，在时间t时刻，作用在每个粒子的力的总和等于它与其他所有粒子的相互作用力的矢量和。根据此力，可以得到此粒子的加速度，结合t时刻的位置与速度，可以得到t+dt时刻的位置与速度。这力在此时间间隔期间假定为为常数。作用在新位置上的粒子的力可以求出，然后可以导出t+2dt时刻的位置与速度等与此类似。分子动力学方程的数值求解有多种方法，各有其优缺点。下面将分别介绍各种常见的计算方法。4.分子运动方程的数值求解Verlet算法Verlet于1967年提出的，在分子动力学中，积分运动方程运用最广泛的方法。这种算法运用t时刻的位置和速度及t-dt时刻的位置，计算出t+dt时刻的位置r(t+dt)：计算速度有很多种方法，一个简单的方法是用t+dt时刻与t-dt时刻的位置差除以2dt：Verlet算法简单，存储要求适度，但它的一个缺点是位置r(t+dt)要通过小项与非常大的两项2r(t)与r(t-dt)的差得到，这容易造成精度损失。4.分子运动方程的数值求解Leap-frog算法Hockney在1970年提出Leap-frog算法。速度表示为：原子坐标的微分可以这样表述为由t-0.5dt时刻的速度与t时刻的加速度计算出速度v(t+0.5dt)，然后计算出位置r(t+dt）。Leap-frog算法相比Ver-let算法有两个优点:它包括显速度项，并且计算量稍小。它也有明显的缺陷：位置与速度不是同步的。4.分子运动方程的数值求解速度Verlet算法Swope在1982年提出的速度Verlet算法可以同时给出位置、速度与加速度，且不牺牲精度：该算法需要储存每个时间步的坐标、速度和力。该方法的每个时间步涉及两个时间步，需要计算坐标更新以后和速度更新前的力。在这三种算法中，速度Verlet算法精度和稳定性最好。力的计算是很费机时,如果内存不是问题，高阶方法可以使用较大的时间步长，因此对计算是非常有利的。4.分子运动方程的数值求解预测-校正算法（Gear算法）为了使用尽可能大的时间步长，或者在相同的时间步长时获得较高的精度，可以储存和使用前一步的力，并使用预测-校正算法更新位置和速度。Gear于1971年提出了基于预测-校正积分方法的算法。首先，根据Taylor展开式，预测新的位置、速度与加速度。然后，根据新预测的位置，计算t+dt时刻的力，然后计算加速度。最后，这加速度再与由Taylor级数展开式预测的加速度进行比较。利用两者之差用来校正位置与速度项。Gear的预测-校正算法需要的存储量为3(O+l)N，O是应用的最高阶微分数，N是原子数目。5.边界条件与初值边界条件由于计算机的运算能力有限，模拟系统的粒子数不可能很大。在数值计算中，对于平衡态分子动力学