分子动力学作业讲解

1/23分子动力学(MD)1分子动力学（MD）基础1.1MD分类1.2MD简介1.3MD适用范围2分子动力学运动方程数值求解2.1基础知识2.1.1运动方程2.1.2空间描述2.1.3最小作用量原理2.1.4拉格朗日(Lagrange)方程2.1.5哈密顿(Hamilton)方程2.2粒子运动方程的数值解法2.2.1Verlet算法2.2.2欧拉(Euler)预测—矫正公式2.2.3Gear预测—矫正方法3分子动力学原胞与边界条件3.1分子动力学原胞3.2边界条件3.2.1自由表面边界3.2.2固定边界3.2.3柔性边界3.2.4周期性边界4势函数与分子力场4.1势函数4.1.1两体势4.1.2多体势4.2分子力场4.2.1分子力场函数的构成2/234.2.2常用力场函数和分类5分子动力学模拟的基本步骤5.1设定模拟所采用的模型5.2给定初始条件5.3趋于平衡计算5.4宏观物理量的计算6平衡态分子动力学模拟6.1系综6.2微正则系综的分子动力学模拟6.3正则系综的分子动力学模拟3/231分子动力学（MD）基础1.1MD分类微正则系综（VNE）正则系综（VNP）平衡态MD等温等压系综（NPT）经典MD等焓等压系综（NPH）巨正则系综（ＶＴμ）非平衡态ＭＤ量子MD1.2分子动力学(MD)简介分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法。分子动力学方法为确定性模拟方法，广泛地用于研究经典的多粒子体系的研究中,是按该体系内部的内禀动力学规律来计算并确定位形的转变。分子动力学方法是通过建立一组分子的运动方程，并通过直接对系统中的一个个分子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计计算方法得到多体系统的静态和动态特性,从而得到系统的宏观性质。在分子动力学中，粒子的运动行为是通过经典的Newton运动方程所描述。系统的所有粒子服从经典力学的运动规律，它的动力学方程就是从经典力学的运动方程——拉格朗日(lagrange)方程和哈密顿(Hamilton)方程导出。1.3适用范围原则上，分子动力学方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分子，也可以是其它的微观粒子。实际上，分子动力学模拟方法和随机模拟方法一样都面临着两个基本限制：4/23一个是有限观测时间的限制；另一个是有限系统大小的限制。通常人们感兴趣的是体系在热力学极限下（即粒子数目趋于无穷时）的性质。但是计算机模拟允许的体系大小要比热力学极限小得多，因此可能会出现有限尺寸效应。为了减小有限尺寸效应，人们往往引入周期性、全反射、漫反射等边界条件。当然边界条件的引入显然会影响体系的某些性质。2分子动力学运动方程数值求解2.1基础知识2.1.1运动方程系统的动力学机制决定运动方程的形式。在分子动力学方法处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学。每个分子运动的内禀动力学是用理论力学上的哈密顿量或者拉格朗日量来描述，也可以直接用牛顿运动方程来描述。采用分子动力学方法时，必须对一组分子运动微分方程做数值求解。从计算数学的角度来看，这是个求一个初值问题的微分方程的解。实际上计算数学为了求解这种问题已经发展了许多的算法。但是并不是所有的这些算法都可以用来解决物理问题。2.1.2空间描述在空间描述如何物体的运动，如果其本身的大小可以忽略时，就可以将其看作是粒子（或质点）。粒子描述：空间位置：r速度：v=dr/dt加速度：若一个系统由N个粒子组成，则粒子描述：空间位置：r1,r2,r3…,rN笛卡尔坐标系，粒子有3N个自由度设系统有s个自由度广义坐标：q1,q2,q3,…,qN5/23广义速度：q1,q2,q3,…,qN2.1.3最小作用量原理莫培督1744年提出最小作用量原理：保守的、完整的力学系统，由某一初位形转变到另一位形的一切具有相同能量的可能运动中，真实的运动是其作用量具有极小值的那种运动。力学系统中，构造能量函数L及其作用量S作用量的积分式叫做泛函(functional)，作用量取极值的方法就是求其变分δS=0。2.1.4拉格朗日(Lagrange)方程由最小作用量原理可导出拉格朗日方程对于孤立的保守系统，每个粒子在势场U中运动，则系统整体的Lagrange函数是得到第i个粒子的牛顿运动方程（α指每个粒子的自由度）2.1.5哈密顿(Hamilton)方程哈密顿(Hamilton)原理:保守的、完整的力学系统在相同时间内，由某一初位形转移到另一已知位形的一切可能运动中，真实运动的作用函数具有极值，即作用函数的变分等于零。6/23哈密顿(Hamilton)方程，Lagrange函数全微分形式：则定义哈密顿函数或哈密顿量为：哈密顿函数H是动量和坐标的函数，是动能和势能之和：变量为动量p和坐标r的Hamilton方程：这就是变量为动量p和坐标q的哈密顿方程。如果系统的哈密顿函数不显含时间，就有dH/dt=0，即得到能量守恒定律。2.2粒子运动方程的数值解法设粒子的坐标、速度、动量及其作用力分别用x(t)，v(t)，p(t)，f(x,t)表示，其初始值为x(0)，v(0)，p(0)，f(0)。则决定粒子运动的牛顿方程是7/23同时由于数值计算求解微分方程是用差分的方法，习惯上将时间的变化间隔Δt用h表示，叫做时间步长。坐标预测的Taylor展开式：势函数U应看作是离子间相互作用势V和外势Uext2.2.1Verlet算法根据Taylor公式，在t时刻求t+h时刻的坐标和作用力时，分成向前和向后的Taylor展开式将上面两式相加得到从而得到坐标的计算公式将两个Taylor展开式相减可以得到8/23从而得到速度或动量的计算公式和加速度或作用力的计算公式于是得到用上一个时间点(t-h)和当前时间点(t)的坐标r和速度v及作用力f，来计算出下一个时间点(t+h)的坐标和速度及作用力的三点公式。这种方法称为Verlet方法。2.2.2欧拉(Euler)预测—矫正公式第一种方法，开放式或单预测式。对下一步的位置和力的计算仅与前几步的已知位置和力有关。实际计算证明，这种计算方式的精度小而误差大，所以实际计算中除了在特殊情况下，已不再采用。第二种方法，封闭式或预测——矫正式。它的计算步骤是：用预测公式计算得到下一步的数值yn+1；将其带入函数式，计算f(yn+1)去矫正预测值yn+1；从而得到矫正后的rn+1。具体操作看下面的欧拉(Euler)预测——矫正公式：预测值矫正值9/232.2.3Gear预测—矫正方法Gear发展出预测－矫正方法(Predict-corrector)。经证明，这是一种精度很高的完全适用于分子动力学的算法被广泛应用。为方便，使用矢量记法。将下一步预测值的每一项进行Taylor展开组成一个列矢量，为书写方便记做转置形式称为N-表象矢量。若将当前和以前几个时刻的坐标值作为一个矢量的各元素，则有称为C-表象矢量。还有一种使用起来比较方便的F-表象，矢量的元素由当前的坐标、速度和当前几前两步的力（或加速度）构成称为F-表象矢量。Gear指出，各种表象可以通过一个变换进行相互转换。矢量rn+1也就是所有n-1阶多项式空间中的矢量，将一个表象R1变换到另一个表象R2的变换完全是一个在该空间中的正交线形变换其中，变换矩阵T可以很容易找到。(1)预测算法。用矩阵代数的形式描述预测值为：yn+1是rn+1的预测值。预测矩阵A是从某一物理量的对时间步长的Taylor展开中得到。N-表象中，5值预测矩阵为：10/23A中的各列元素值构成从0阶到5阶的二项式的展开系数。更高阶的预测矩阵只需在矩阵最后一列添加更高阶二项式的展开系数。预测矩阵的变换矩阵为由此不难得到预测矩阵在其他表象中的形式。一个4值F-表象中的预测矩阵为(2)矫正算法。为了得到r的下一步预测值rk+1，需要将矢量rk+1的每一项进行Taylor展开，由于t=kh，得到rk+1=Brk，B是(n+1)×(n+1)的系数矩阵，其第i行第j列的元素值为把预测项和矫正项结合起来C是矫正系数矢量。通过求解上面方程的本征值，即可求出系数向量C的各元素值。如果Taylor展开的项数不太多，C的各值也可以直接用代入法求得。Gear推导的4~8值N-表象中矫正矢量C的值为11/23代人预测值yn+1，计算预测位置上的加速度f(yn+1)/m，得到更为准确的下一步位置可以由矫正式：3分子动力学原胞与边界条件3.1分子动力学原胞分子动力学模拟方法往往用于研究大块物质在给定密度下的性质，而实际计算模拟不可能在几乎是无穷大的系统中进行。所以必须引进一个叫做分子动力学原胞的体积元，以维持一个恒定的密度。对气体和液体，如果所占体积足够大，并且系统处于热平衡状态的情况下，那么这个体积的形状是无关紧要的。对于晶态的系统，原胞的形状是有影响的。为了计算简便，取一个立方形的体积为分子动力学原胞。设原胞的线度大小为L，则体积为L3。由于引进这样的立方体箱子，将产生六个我们不希望出现的表面。模拟中碰撞这些箱子的表面的粒子应当被反射回到原胞内部，特别是对粒子数目很少的系统。然而这些表面的存在对系统的任何一种性质都会有重大的影响。12/233.2边界条件采用分子动力学方法，必需对被计算的粒子系统给定适当的边界条件。这些边界条件大致可分成四种。（1）自由表面边界(free-surfaceboundary)这种边界条件常用于大型的自由分子的模拟。（2）固定边界(rigidboundary)在所有要计算到粒子晶胞之外还要包上几层结构相同的位置不变的粒子，包层的厚度必须大于粒子间相互作用的力程范围。包层部分代表了与运动粒子起作用的宏观晶体的那一部分。（3）柔性边界(flexibleboundary)这种边界比固定边界更接近实际。它允许边界上的粒子有微小的移动以反映内层粒子的作用力施加到它们身上时的情况。（4）周期性边界(periodicboundary)13/23在模拟较大的系统时，为了消除表面效应或边界效应，常采用周期性边界条件。就是让原胞上下、左右、前后对边上的粒子间有相互作用。为了将分子动力学原胞有限立方体内的模拟，扩展到真实大系统的模拟，通常采用周期性边界条件。采用周期性边界条件，就可以消除引入原胞后的表面效应，构造出一个准无穷大的体积来更精确地代表宏观系统。实际上，这里我们做了一个假定，即让这个小体积原胞镶嵌在一个无穷大的大块物质之中。做一个假设：让这个小体积原胞镶嵌在一个无穷大的大块物质中，周期性边界条件的数学表示形式为：其中A为任意可观测量，n1,n2,n3为任意整数，这个边界条件就是命令基本MD原胞完全等同地重复无穷多次。4势函数与分子力场4.1势函数由于物质系统的复杂性以及原子间相互作用类型的不同，很难得到满足各种不同体系和物质的一般性而又精度很高的势函数，所以针对不同的物质体系人们陆续发展了大量的经验和半经验的势函数。半经验势：形式经过一定的理论分析而得到，但其中一些参数则需要根据宏观实验参数用经验方法来确定，这些宏观实验参数主要有弹性常数、平衡点阵常数以及内聚能、空位形成能和层错能等。经验势：为了拟合的方便，在选择势函数的形式时，并不一定要求有确切的理论依据，而是出于经验的估计和拟合方便的需要，相对自由的选择势函数的形式。14/23势函数一般分为二体势能和三体以上的多体势能函数。势函数形式：4.1.1两体势两体势认为原子之间的相互作用是两两原子之间的作用，与其他原子的位置无关。这种势的形式可以比较充分地描述除半导体、金属以外的无机化合物材料中的相互作用。（1）硬球模型（2）Lennard-Jones(LJ)势（3）Morse势15/23（4）Borna-Mayer势4.1.2多体势两体势模型在解决分子晶体和离子型化合物问题中取得了成功，但对于含有共价键的过渡金属却遇到了很多困难。（1）嵌入原子势（EAM，Em