分子动力学方法_zkw

分子动力学方法MolecularDynamicsSimulations张凯旺2009年02月18日湘潭大学物理学专业课材料与光电物理学院目录第一节引言第二节分子动力学运动方程数值求解第三节分子动力学原胞与边界条件第四节势函数与分子力场第五节分子动力学模拟的基本步骤第六节平衡态分子动力学模拟第一节引言•分子动力学方法（MolecularDynamics，简称MD）Alder和Wainght在1957年至1959年间应用于理想“硬球”液体模型，很多简单液体中分子之间的相互作用的重要性质在两人的研究中被发现；Rahman于1964年应用一种更接近的液体模型模拟了液氦；Verlet于1967~1971年，Levesque于1983年做了进一步的研究，模拟了更复杂的液体——水、溶盐、聚合物和蛋白质等；使用分子动力学模拟研究真实系统在1974年由Rahman和Stilliger完成，研究对象为液体水；第一个蛋白质的分子动力学研究报道是McCammon在1977年，研究对象是牛胰岛素抑制剂。根本问题•分子动力学计算机模拟是在原子量级上模拟材料的性质，模拟的根本问题是要确定一群有相互作用的粒子在时空中的演化规律，也就是说，要知道各个粒子什么时候在什么地方？是如何运动的？首先要建立数学模型，即把关于微观粒子或粒子团的结构、粒子间力的知识与牛顿力学结合起来，指定粒子运动应遵循的自然规律和粒子间的相互作用的形式；然后用计算机计算粒子集合的相轨迹，从而确定系统的静态和动态性质。计算机模拟分类•对于一个多粒子体系的实验观测物理量的数值可以由总的平均得到。•但是由于实验体系又非常大，不可能计算求得所有涉及到的态的物理量数值的总平均。•按照产生位形变化的方法，有两类方法对有限的一系列态的物理量做统计平均。随机模拟方法分子动力学方法计算机模拟分类•随机模拟方法。随机模拟方法是实现Gibbs的统计力学途径。在随机模拟方法中，体系位形的转变是通过马尔科夫(Markov)过程，由随机性的演化引起的。这里的马尔科夫过程相当于是内禀动力学在概率方面的对应物。该方法可以被用到没有任何内禀动力学模型体系的模拟上。优点:随机模拟方法计算的程序简单，占内存少，但是该方法难于处理非平衡态的问题。计算机模拟分类•分子动力学方法。分子动力学方法为确定性模拟方法，广泛地用于研究经典的多粒子体系的研究中。分子动力学方法是按该体系内部的内禀动力学规律来计算并确定位形的转变。首先需要建立一组分子的运动方程，并通过直接对系统中的一个个分子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计计算方法得到多体系统的静态和动态特性,从而得到系统的宏观性质。计算机模拟分类分子动力学模拟方法可以看作是体系在一段时间内的发展过程的模拟。在这样的处理过程中可以看出：分子动力学方法中不存在任何随机因素。在分子动力学方法处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的。系统的动力学机制决定运动方程的形式。计算机模拟分类在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学。每个分子运动的内禀动力学是用理论力学上的哈密顿量或者拉格朗日量来描述，也可以直接用牛顿运动方程来描述。这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。但是使用该方法的程序较复杂，计算量大，占内存也多。分子动力学简介•分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法。•通过求解所有粒子的运动方程，分子动力学方法可以用于模拟与粒子运动路径相关的基本过程。•在分子动力学中，粒子的运动行为是通过经典的Newton运动方程所描述。系统的所有粒子服从经典力学的运动规律，它的动力学方程就是从经典力学的运动方程——拉格朗日(lagrange)方程和哈密顿(Hamilton)方程导出。分子动力学简介•分子动力学得天独厚的优势是能够模拟分子的运动轨道，即能够提供分子运动和变化的最为微观的、细致的信息。•这对传统的力学、热力学、光谱学等实验方法都是难以办到的。•即使是现代分析理论和实验测量手段，要想准确的描述物质的非平衡态性质都是非常困难的。分子动力学简介•但是如果有分子动力学计算得到的运行轨道，能通过非常直接的计算公式给出感兴趣的非平衡态物理量。•这种特点对现今所有的理论、实验、计算手段来讲都是独一无二的。•正是这种不可替代的作用，使得分子动力学在当今自然科学研究中获得了如此广泛和深入的应用。适用范围广泛•原则上，分子动力学方法所适用的微观物理体系并无什么限制。•这个方法适用的体系既可以是少体系统，也可以是多体系统；•既可以是点粒子体系，也可以是具有内部结构的体系；•处理的微观客体既可以是分子，也可以是其它的微观粒子。适用范围广泛•自五十年代中期开始，分子动力学方法得到了广泛的应用。•它与蒙特卡洛方法一起已经成为计算机模拟的重要方法。应用分子动力学方法取得了许多重要成果。•例如气体或液体的状态方程、相变问题、吸附问题等，以及非平衡过程的研究。•其应用已从化学反应、生物学的蛋白质到重离子碰撞等广泛的学科研究领域。实际使用的限制•实际上，分子动力学模拟方法和随机模拟方法一样都面临着两个基本限制：一个是有限观测时间的限制；另一个是有限系统大小的限制。通常人们感兴趣的是体系在热力学极限下（即粒子数目趋于无穷时）的性质。但是计算机模拟允许的体系大小要比热力学极限小得多，因此可能会出现有限尺寸效应。为了减小有限尺寸效应，人们往往引入周期性、全反射、漫反射等边界条件。当然边界条件的引入显然会影响体系的某些性质。实际使用的限制•对体系的分子运动方程组采用计算机进行数值求解时，需要将运动方程离散化为有限差分方程。常用的求解方法有欧拉法、龙格-库塔法、辛普生法等。数值计算的误差阶数显然取决于所采用的数值求解方法的近似阶数。原则上，只要计算机计算速度足够大，内存足够多，我们可以使计算误差足够小。第二节分子动力学运动方程数值求解•系统的动力学机制决定运动方程的形式。•在分子动力学方法处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的。•在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学。•每个分子运动的内禀动力学是用理论力学上的哈密顿量或者拉格朗日量来描述，也可以直接用牛顿运动方程来描述。•这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。运动方程•采用分子动力学方法时，必须对一组分子运动微分方程做数值求解。•从计算数学的角度来看，这是个求一个初值问题的微分方程的解。•实际上计算数学为了求解这种问题已经发展了许多的算法。•但是并不是所有的这些算法都可以用来解决物理问题。空间描述•在空间描述如何物体的运动，如果其本身的大小可以忽略时，就可以将其看作是粒子（或质点）。•粒子描述：空间位置：r速度：v=dr/dt加速度：a=dv/dt=d2r/dt2•若一个系统由N个粒子组成，则粒子描述：空间位置：r1，r2，…，rN笛卡尔坐标系，粒子有3N个自由度•设系统有s个自由度广义坐标：q1，q2，…，qs广义速度：sqqq,,,21最小作用量原理•莫培督1744年提出最小作用量原理：保守的、完整的力学系统，由某一初位形转变到另一位形的一切具有相同能量的可能运动中，真实的运动是其作用量具有极小值的那种运动。•力学系统中，构造能量函数L及其作用量S•作用量的积分式叫做泛函(functional)，作用量取极值的方法就是求其变分δS=0。0d,,d,,,,2121tttttttqtqLStttqtqLtqSttqtqLL拉格朗日(Lagrange)方程•由最小作用量原理可导出拉格朗日方程•对于孤立的保守系统，每个粒子在势场U中运动，则•系统整体的Lagrange函数是•得到第i个粒子的牛顿运动方程（α指每个粒子的自由度）iiiiiiiiiiiiiiiiiiiiUqUqmUqmtqqLLUqmtqqLqLqLt-----aa2221,,21,,0dd哈密顿(Hamilton)方程•哈密顿(Hamilton)原理:保守的、完整的力学系统在相同时间内，由某一初位形转移到另一已知位形的一切可能运动中，真实运动的作用函数具有极值，即作用函数的变分等于零。•哈密顿(Hamilton)方程：Lagrange函数全微分形式：因为则有哈密顿(Hamilton)方程由于有整理得定义哈密顿函数或哈密顿量哈密顿(Hamilton)方程哈密顿函数H是动量和坐标的函数，是动能和势能之和变量为动量p和坐标r的Hamilton方程这就是变量为动量p和坐标q的哈密顿方程。如果系统的哈密顿函数不显含时间，就有dH/dt=0，即得到能量守恒定律。粒子运动方程的数值解法设粒子的坐标、速度、动量及其作用力分别用x(t)，v(t)，p(t)，f(x,t)表示，其初始值为x(0)，v(0)，p(0)，f(0)。则决定粒子运动的牛顿方程是同时粒子运动方程的数值解法由于数值计算求解微分方程是用差分的方法，习惯上将时间的变化间隔∆t用h表示，叫做时间步长。坐标预测的Taylor展开式：势函数U应看作是离子间相互作用势V和外势Uext之和Verlet算法根据Taylor公式，在t时刻求t+h时刻的坐标和作用力时，分成向前和向后的Taylor展开式将上面两式相加得到从而得到坐标的计算公式Verlet算法将两个Taylor展开式相减可以得到从而得到速度或动量的计算公式和加速度或作用力的计算公式于是得到用上一个时间点(t-h)和当前时间点(t)的坐标r和速度v及作用力f，来计算出下一个时间点(t+h)的坐标和速度及作用力的三点公式。这种方法称为Verlet方法。欧拉(Euler)预测—矫正公式•计算方法至少应当符合对力的一阶导数的恰当处理，以免误差的积累和解的不稳定。•在相关坐标系中求解方程时，以时间步长h的阶数不应小于3。•采用更高的阶数则要求作用力的更高阶导数，这就有可能带来不确定性和非系统性的问题。•复杂系统中的计算常带有截断误差，使得高于3阶或4阶的算法并不具有更大的优越性，因阶数越高，数值会越小，在与数值相对较大的低阶项相加减时会被大数“吃掉”。欧拉(Euler)预测—矫正公式第一种方法，开放式或单预测式。对下一步的位置和力的计算仅与前几步的已知位置和力有关。实际计算证明，这种计算方式的精度小而误差大，所以实际计算中除了在特殊情况下，已不再采用。第二种方法，封闭式或预测——矫正式。它的计算步骤是：用预测公式计算得到下一步的数值yn+1；将其带入函数式，计算f(yn+1)去矫正预测值yn+1；从而得到矫正后的rn+1。欧拉(Euler)预测—矫正公式•具体操作看下面的欧拉(Euler)预测——矫正公式：预测值矫正值Gear预测—矫正方法•Gear发展出预测－矫正方法(Predict-corrector)。经证明，这是一种精度很高的完全适用于分子动力学的算法，被广泛应用。•为方便，使用矢量记法。将下一步预测值的每一项进行Taylor展开•组成一个列矢量，为书写方便记做转置形式称为N-表象矢量Gear预测—矫正方法•若将当前和以前几个时刻的坐标值作为一个矢量的各元素，则有称为C-表象矢量•还有一种使用起来比较方便的F-表象，矢量的元素由当前的坐标、速度和当前几前两步的力（或加速度）构成称为F-表象矢量Gear预测—矫正方法•Gear指出，各种表象可以通过一个变换进行相互转换。•矢量rn+1也就是所有n-1阶多项式空间中的矢量，将一个表象R1变换到另一个表象R2的变换完全是一个在该空间中的正交线形变换其中，变换矩阵T可以很容易找到。Gear预测—矫正方法•比如一个3值的N-表象中的矢量变换到C-表象中。由Taylor展开得到则就有表象转换运算式Gear预测—矫正方法•同样，一个4值F-表象的矢量可以从4值N-表象变换得到G