分子动力学(MD) 1 分子动力学（MD）基础
1.1 MD分类
1.2 MD简介
1.3 MD适用范围
2 分子动力学运动方程数值求解
2.1 基础知识
2.1.1 运动方程
2.1.2 空间描述
2.1.3 最小作用量原理
2.1.4 拉格朗日(Lagrange)方程
2.1.5 哈密顿(Hamilton)方程
2.2 粒子运动方程的数值解法
2.2.1 Verlet算法
2.2.2 欧拉(Euler)预测—矫正公式
2.2.3 Gear预测—矫正方法
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 分子力场函数的构成
4.2.2 常用力场函数和分类
5 分子动力学模拟的基本步骤
5.1 设定模拟所采用的模型
5.2 给定初始条件
5.3 趋于平衡计算
5.4 宏观物理量的计算
6 平衡态分子动力学模拟
6.1 系综
6.2 微正则系综的分子动力学模拟6.3 正则系综的分子动力学模拟
1 分子动力学（MD）基础
1.1MD分类
微正则系综（VNE）
正则系综（VNP）
平衡态MD 等温等压系综（NPT）
经典MD 等焓等压系综（NPH）
巨正则系综（ＶＴμ）
非平衡态ＭＤ
量子MD
1.2分子动力学(MD)简介
分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法。分子动力学方法为确定性模拟方法，广泛地用于研究经典的多粒子体系的研究中,是按该体系内部的内禀动力学规律来计算并确定位形的转变。
分子动力学方法是通过建立一组分子的运动方程，并通过直接对系统中的一个个分子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计计算方法得到多体系统的静态和动态特性, 从而得到系统的宏观性质。
在分子动力学中，粒子的运动行为是通过经典的Newton运动方程所描述。系统的所有粒子服从经典力学的运动规律，它的动力学方程就是从经典力学的运动方程——拉格朗日(lagrange)方程和哈密顿(Hamilton)方程导出。
1.3适用范围
原则上，分子动力学方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分子，也可以是其它的微观粒子。
实际上，分子动力学模拟方法和随机模拟方法一样都面临着两个基本限制：
一个是有限观测时间的限制；另一个是有限系统大小的限制。通常人们感兴趣的是体系在热力学极限下（即粒子数目趋于无穷时）的性质。但是计算机模拟允许的体系大小要比热力学极限小得多，因此可能会出现有限尺寸效应。为了减小有限尺寸效应，人们往往引入周期性、全反射、漫反射等边界条件。当然边界条件的引入显然会影响体系的某些性质。
2 分子动力学运动方程数值求解
2.1 基础知识
2.1.1 运动方程
系统的动力学机制决定运动方程的形式。在分子动力学方法处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学。每个分子运动的内禀动力学是用理论力学上的哈密顿量或者拉格朗日量来描述，也可以直接用牛顿运动方程来描述。
采用分子动力学方法时，必须对一组分子运动微分方程做数值求解。从计算数学的角度来看，这是个求一个初值问题的微分方程的解。实际上计算数学为了求解这种问题已经发展了许多的算法。但是并不是所有的这些算法都可以用来解决物理问题。
2.1.2 空间描述
在空间描述如何物体的运动，如果其本身的大小可以忽略时，就可以将其看作是粒子（或质点）。
粒子描述：空间位置：r
速度：v = dr/dt
加速度：
若一个系统由N个粒子组成，则粒子描述：
空间位置：r1,r2,r3…,r N
笛卡尔坐标系，粒子有3N个自由度
设系统有s个自由度
广义坐标：q1,q2,q3,…,q N
广义速度：q1,q2,q3,…,q N
2.1.3 最小作用量原理
莫培督1744年提出最小作用量原理：保守的、完整的力学系统，由某一初位形转变到另一位形的一切具有相同能量的可能运动中，真实的运动是其作用量具有极小值的那种运动。
力学系统中，构造能量函数L及其作用量S
作用量的积分式叫做泛函(functional)，作用量取极值的方法就是求其变分δS = 0。
2.1.4 拉格朗日(Lagrange)方程
由最小作用量原理可导出拉格朗日方程
对于孤立的保守系统，每个粒子在势场U中运动，则
系统整体的Lagrange函数是
得到第i个粒子的牛顿运动方程（α指每个粒子的自由度）
2.1.5 哈密顿(Hamilton)方程
哈密顿(Hamilton)原理:保守的、完整的力学系统在相同时间内，由某一初位形转移到另一已知位形的一切可能运动中，真实运动的作用函数具有极值，即作用函数的变分等于零。
哈密顿(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)。则决定粒子运动的牛顿方程是
同时
由于数值计算求解微分方程是用差分的方法，习惯上将时间的变化间隔Δt 用h 表示，叫做时间步长。
坐标预测的Taylor展开式：
势函数U 应看作是离子间相互作用势V 和外势Uext
2.2.1 Verlet算法
根据Taylor公式，在t 时刻求t +h 时刻的坐标和作用力时，分成向前和向后的Taylor展开式
将上面两式相加得到
从而得到坐标的计算公式
将两个Taylor展开式相减可以得到
从而得到速度或动量的计算公式
和加速度或作用力的计算公式
于是得到用上一个时间点(t-h)和当前时间点(t)的坐标r和速度v及作用力f，来计算出下一个时间点(t+h)的坐标和速度及作用力的三点公式。这种方法称为Verlet方法。
2.2.2 欧拉(Euler)预测—矫正公式
第一种方法，开放式或单预测式。对下一步的位置和力的计算仅与前几步的已知位置和力有关。实际计算证明，这种计算方式的精度小而误差大，所以实际计算中除了在特殊情况下，已不再采用。
第二种方法，封闭式或预测——矫正式。它的计算步骤是：用预测公式计算得到下一步的数值y n+1；将其带入函数式，计算f(y n+1)去矫正预测值y n+1；从而得到矫正后的r n+1。
具体操作看下面的欧拉(Euler)预测——矫正公式：
预测值
矫正值
2.2.3 Gear预测—矫正方法
Gear发展出预测－矫正方法(Predict-corrector)。经证明，这是一种精度很高的完全适用于分子动力学的算法被广泛应用。为方便，使用矢量记法。
将下一步预测值的每一项进行Taylor展开
组成一个列矢量，为书写方便记做转置形式
称为N-表象矢量。
若将当前和以前几个时刻的坐标值作为一个矢量的各元素，则有
称为C-表象矢量。
还有一种使用起来比较方便的F-表象，矢量的元素由当前的坐标、速度和当前几前两步的力（或加速度）构成
称为F-表象矢量。
Gear指出，各种表象可以通过一个变换进行相互转换。矢量rn+1也就是所有n-1阶多项式空间中的矢量，将一个表象R1变换到另一个表象R2的变换完全是一个在该空间中的正交线形变换
其中，变换矩阵T可以很容易找到。
(1) 预测算法。用矩阵代数的形式描述预测值为：
y n+1是r n+1的预测值。
预测矩阵A是从某一物理量的对时间步长的Taylor展开中得到。N-表象中，5值预测矩阵为：
A中的各列元素值构成从0阶到5阶的二项式的展开系数。更高阶的预测矩阵只需在矩阵最后一列添加更高阶二项式的展开系数。
预测矩阵的变换矩阵为
由此不难得到预测矩阵在其他表象中的形式。一个4值F-表象中的预测矩阵为
(2) 矫正算法。
为了得到r 的下一步预测值r k+1，需要将矢量r k+1的每一项进行Taylor展开，由于t=kh，得到r k+1=B r k，B是(n+1)×(n+1)的系数矩阵，其第i 行第j 列的元素值为
把预测项和矫正项结合起来
C是矫正系数矢量。通过求解上面方程的本征值，即可求出系数向量C的各元素值。如果Taylor展开的项数不太多，C的各值也可以直接用代入法求得。
Gear推导的4~8值N-表象中矫正矢量C的值为
代人预测值yn+1，计算预测位置上的加速度f(yn+1)/m，得到更为准确的下一步位置可以由矫正式：
3 分子动力学原胞与边界条件
3.1 分子动力学原胞
分子动力学模拟方法往往用于研究大块物质在给定密度下的性质，而实际计算模拟不可能在几乎是无穷大的系统中进行。所以必须引进一个叫做分子动力学原胞的体积元，以维持一个恒定的密度。
对气体和液体，如果所占体积足够大，并且系统处于热平衡状态的情况下，那么这个体积的形状是无关紧要的。对于晶态的系统，原胞的形状是有影响的。
为了计算简便，取一个立方形的体积为分子动力学原胞。设原胞的线度大小为L，则体积为L3。由于引进这样的立方体箱子，将产生六个我们不希望出现的表面。模拟中碰撞这些箱子的表面的粒子应当被反射回到原胞内部，特别是对粒子数目很少的系统。然而这些表面的存在对系统的任何一种性质都会有重大的影响。
3.2 边界条件
采用分子动力学方法，必需对被计算的粒子系统给定适当的边界条件。这些边界条件大致可分成四种。
（1）自由表面边界(free-surface boundary)
这种边界条件常用于大型的自由分子的模拟。
（2）固定边界(rigid boundary)
在所有要计算到粒子晶胞之外还要包上几层结构相同的位置不变的粒子，包层的厚度必须大于粒子间相互作用的力程范围。包层部分代表了与运动粒子起作用的宏观晶体的那一部分。
（3）柔性边界(flexible boundary)
这种边界比固定边界更接近实际。它允许边界上的粒子有微小的移动以反映内层粒子的作用力施加到它们身上时的情况。
（4）周期性边界(periodic boundary)
在模拟较大的系统时，为了消除表面效应或边界效应，常采用周期性
边界条件。就是让原胞上下、左右、前后对边上的粒子间有相互作用。
为了将分子动力学原胞有限立方体内的模拟，扩展到真实大系统的模拟，通常采用周期性边界条件。
采用周期性边界条件，就可以消除引入原胞后的表面效应，构造出一个准无穷大的体积来更精确地代表宏观系统。实际上，这里我们做了一个假定，即让这个小体积原胞镶嵌在一个无穷大的大块物质之中。
做一个假设：让这个小体积原胞镶嵌在一个无穷大的大块物质中，周期性边界条件的数学表示形式为：
其中A为任意可观测量，n1,n2,n3为任意整数，这个边界条件就是命令基本MD原胞完全等同地重复无穷多次。
4 势函数与分子力场
4.1 势函数
由于物质系统的复杂性以及原子间相互作用类型的不同，很难得到满足各种不同体系和物质的一般性而又精度很高的势函数，所以针对不同的物质体系人们陆续发展了大量的经验和半经验的势函数。
半经验势：形式经过一定的理论分析而得到，但其中一些参数则需要根据宏观实验参数用经验方法来确定，这些宏观实验参数主要有弹性常数、平衡点阵常数以及内聚能、空位形成能和层错能等。
经验势：为了拟合的方便，在选择势函数的形式时，并不一定要求有确切的理论依据，而是出于经验的估计和拟合方便的需要，相对自由的选择势函数的形式。
势函数一般分为二体势能和三体以上的多体势能函数。
势函数形式：
4.1.1 两体势
两体势认为原子之间的相互作用是两两原子之间的作用，与其他原子的位置无关。这种势的形式可以比较充分地描述除半导体、金属以外的无机化合物材料中的相互作用。
（1）硬球模型
（2）Lennard-Jones(LJ)势
（3）Morse势
（4）Borna-Mayer势
4.1.2 多体势
两体势模型在解决分子晶体和离子型化合物问题中取得了成功，但对于含有共价键的过渡金属却遇到了很多困难。
（1）嵌入原子势（EAM，Embedded atom model）
（2）修正的嵌入原子势（MEAM势）
以上讨论的势都是中心对称的，对于由共价键结合的有机分子和半导体材料是不适宜的。为了更好的描述各种含有共价键作用的材料，发展了很多考虑角度效应的多体势，其中修正的嵌入原子势就是比较常用的一种，它的形式与嵌入原子势一样，只是在背景电子密度中考虑了角度效应，其形式为：
（3）Stillinger-Weber势
4.2分子力场
根据量子力学的波恩-奥本海默近似，一个分子的能量可以近似看作构成分子的各个原子的空间坐标的函数。简单地讲就是分子的能量随分子构型的变化而变化，描述这种分子能量和分子结构之间关系的就是分子力场函数。
分子力场函数为来自实验结果的经验公式，可以讲对分子能量的模拟比较粗糙，但是相比于精确的量子力学从头计算方法，分子力场方法的计算量要小数十倍。而且在适当的范围内，分子力场方法的计算精度与量子化学计算相差无几，因此对大分子复杂体系而言，分子力场方法是一套行之有效的方法。以分子力场
为基础的分子力学计算方法在分子动力学、蒙特卡罗、分子对接等分子模拟方法中有着广泛的应用。
4.2.1 分子力场函数的构成
一般而言，分子力场函数由以下几个部分构成：
键伸缩能：构成分子的各个化学键在键轴方向上的伸缩运动所引起的能量变化。键角弯曲能：键角变化引起的分子能量变化。
二面角扭曲能：单键旋转引起分子骨架扭曲所产生的能量变化。
非键相互作用：包括范德华力、静电相互作用等与能量有关的非键相互作用。交叉能量项：上述作用之间耦合引起的能量变化。
构成一套力场函数体系需要有一套联系分子能量和构型的函数；还需要给出各种不同原子在不同成键状况下的物理参数，比如正常的键长、键角、二面角等；这些力场参数多来自实验或者量子化学计算。
势函数形式, 但一般总表达为分子内与分子间势能之和:
分子内势能(键合)包括键伸缩、键角弯曲和二面角扭转势能；分子间势能(非键合)包括范德华势和静电势, 有的还包括H键：
键合势函数中，一些力场还包含交叉项，使精度更高。交叉项的含义：如键长变化时，键角弯曲势能随键长的不同而不同。
4.2.2 常用力场函数和分类
不同的分子力场会选取不同的函数形式来描述能量与体系构型之间的关系。
到目前，不同的科研团队设计了很多适用于不同体系的力场函数，根据他们选择的函数和力场参数，可以分为以下几类：
传统力场
第二代力场----第二代的势能函数形式比传统力场要更加复杂，涉及的力场参数更多，计算量也更大，当然也相应地更加准确。
通用力场----通用力场也叫基于规则的力场，它所应用的力场参数是基于原子性质计算所得，用户可以通过自主设定一系列分子作为训练集来生成合用的力场参数。
5 分子动力学模拟的基本步骤
5.1 设定模拟所采用的模型
模型的设定，也就是势函数的选取。势函数的研究和物理系统上对物质的描述研究息息相关。
设定模型是分子动力学模拟的第一步工作。例如在一个分子系统中，假定两个分子间的相互作用势为硬球势，其势函数表示为
常用的是Lennard-Jones 势函数、EAM势函数等，不同的物质状态描述用不同的势函数。
5.2 给定初始条件
运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。一般讲，系统的初始条件不可能知道，实际上也不需要精确选择代求系统的初始条件，因为模拟实践足够长时，系统就会忘掉初始条件。当然，合理的初始条件可以加快系统趋于平衡的时间和步伐，获得好的精度。
常用的初始条件可以选择为：令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到；令初始位置随机的偏离差分划分网格的格子上，初始速度为零；令初始位置随机的偏离差分划分网格的格子上，初始速度也是从玻尔兹曼分布随机抽样得到。
5.3 趋于平衡计算
按照给出的运动方程、边界条件和初始条件，就可以进行分子动力学模拟计算。但是，这样计算出的系统不会具有所要求的系统能量，并且这个状态本身也还不是一个平衡态。
为了使系统达到平衡，模拟中需要一个趋衡过程。在这个过程中，增加或从系统中移出能量，直到系统具有所要求的能量。然后，再对运动方程中的时间向前积分若干步，使系统持续给出确定能量值。称这时系统已经达到平衡态。这段达到平衡所需的时间称为弛豫时间。
分子动力学计算的基本思想是：赋予分子体系初始运动状态之后，利用分子的自然运动在相空间中抽取样本进行统计分子动力学计算，时间步长就是抽样的间隔，因而时间步长的选取对动力学模拟非常重要。它决定了模拟所需要的时间。
时间步长h太长会造成分子间激烈碰撞，体系数据溢出；为了减小误差，步长h 必须取得小一些；但是取得太小，系统模拟的弛豫时间就很长，太短的时间步长还会降低模拟过程搜索相空间的能力。因此一般选取的时间步长为体系各个自由度中最短运动周期的十分之一。这里需要积累一定的模拟经验，选择适当的时间步长h 。
5.4 宏观物理量的计算
实际计算宏观物理量往往是在分子动力学模拟的最后阶段进行的。它是沿着相空间轨迹求平均来计算得到的。
如，在模拟过程中计算出的动能值是在不连续的路径上的值，在时间的各个间断点μ上计算动能的平均值
在分子动力学模拟过程中，温度是需要加以监测的物理量，特别是在模拟的起始阶段。根据能量均分定理，我们可以从平均动能值计算得到温度值，
其中d为每个粒子的自由度，如果不考虑系统所受的约束，则d＝3。
6 平衡态分子动力学模拟
在经典分子动力学模拟方法的应用当中，存在着对两种系统状态的分子动力学模拟。一种是对平衡态的MD模拟；另一类是对非平衡态的分子动力学模拟。
对平衡态系综分子动力学模拟又可以分为如下类型：
微正则系综的分子动力学(N,V,E）模拟
正则系综的分子动力学（N,V,T）模拟
等温等压系综分子动力学(N,P,T) 模拟
等焓等压系综分子动力学(N,P,H) 模拟
巨正则系综分子动力学(V,T, μ) 模拟
6.1 系综
系综(ensemble) 是指在一定的宏观条件下(约束条件)，大量性质和结构完全相同的、处于各种运动状态的、各自独立的系统的集合。全称为统计系综。系综是用统计方法描述热力学系统的统计规律性时引入的一个基本概念；系综是统计理论的一种表述方式，系综理论使统计物理成为普遍的微观统计理论；系综并不是实际的物体，构成系综的系统才是实际物体。约束条件是由一组外加宏观参量来表示。在平衡统计力学范畴下，可以用来处理稳定系综。
根据宏观约束条件，系综被分为以下几种：
（1）正则系综(canonical ensemble)，全称应为“宏观正则系综”，简写为NVT，即表示具有确定的粒子数(N)、体积(V)、温度(T)。
正则系综是蒙特卡罗方法模拟处理的典型代表。假定N个粒子处在体积为V 的盒子内，将其埋入温度恒为T的热浴中。此时，总能量(E)和系统压强(P)可能在某一平均值附近起伏变化。
平衡体系代表封闭系统，与大热源热接触平衡的恒温系统。正则系综的特征函数是亥姆霍兹自由能F(N,V,T)。
（2）微正则系综(micro-canonical ensemble)，简写为NVE，即表示具有确定的粒子数(N)、体积(V)、总能量(E)。
微正则系综广泛被应用在分子动力学模拟中。假定N个粒子处在体积为V 的盒子内，并固定总能量(E)。此时，系综的温度(T)和系统压强(P)可能在某一平
《装备制造技术》 2007年第 10期 收稿日期 :2007-08-21 作者简介 :申海兰 , 24岁 , 女 , 河北人 , 在读研究生 , 研究方向为微机电系统。 分子动力学模拟方法概述 申海兰 , 赵靖松 (西安电子科技大学机电工程学院 , 陕西西安 710071 摘要 :介绍了分子动力学模拟的基本原理及常用的原子间相互作用势 , 如Lennard-Jones 势 ; 论述了几种常用的有限差分算法 , 如 Verlet 算法 ; 说明了分子动力学模拟的几种系综及感兴趣的宏观统计量的提取。关键词 :分子动力学模拟 ; 原子间相互作用势 ; 有限差分算法 ; 系综中图分类号 :O3 文献标识码 :A 文章编号 :1672-545X(200710-0029-02 从统计物理学中衍生出来的分子动力学模拟方法 (molec- ular dynamics simulation , M DS , 实践证明是一种描述纳米科技 研究对象的有效方法 , 得到越来越广泛的重视。所谓分子动力学模拟 , 是指对于原子核和电子所构成的多体系统 , 用计算机模拟原子核的运动过程 , 从而计算系统的结构和性质 , 其中每一个原子核被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动 [1]。它被认为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段 , 称之为“计算机实验” 手段 [2], 在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。
根据模拟对象的不同 , 将它分为平衡态分子动力学模拟 (EM DS (和非平衡态分子动力学模拟 (NEM DS 。其中 , EM DS 是分子动力学模拟的基础 ; NEM DS 适用于非线性响应系统的模拟 [3]。下面主要介绍 EM DS 。 1分子动力学方法的基本原理 计算中根据以下基本假设 [4]: (1 所有粒子的运动都遵循经典牛顿力学规律。 (2 粒子之间的相互作用满足叠加原理。 显然这两条忽略了量子效应和多体作用 , 与真实物理系统存在一定差别 , 仍然属于近似计算。 假设 N 为模拟系统的原子数 , 第 i 个原子的质量为 m i , 位置坐标向量为 r i , 速度为 v i =r ? i , 加速度为 a i =r ?? i , 受到的作用力为 F i , 原子 i 与原子 j 之间距离为 r ij =r i -r j , 原子 j 对原子 i 的作用力为 f ij , 原子 i 和原子 j 相互作用势能为 ! (r ij , 系统总的势能为 U (r 1, r 2, K r N = N i =1! j ≠ i ! " (r ij , 所有的物理量都是随时 间变化的 , 即 A=A (t , 控制方程如下 : m i r ?? i =F i =j ≠ i
分子动力学的模拟过程 分子动力学模拟作为一种应用广泛的模拟计算方法有其自身特定的模拟步骤,程序流程也相对固定。本节主要就分子动力学的模拟步骤和计算程序流程做一些简单介绍。 1. 分子动力学模拟步驟 分子动力学模拟是一种在微观尺度上进行的数值模拟方法。这种方法既可以得到一些使用传统方法,热力学分析法等无法获得的微观信息,又能够将实际实验研究中遇到的不利影响因素回避掉,从而达到实验研宄难以实现的控制条件。 分子动力学模拟的步骤为: (1)选取所要研究的系统并建立适当的模拟模型。 (2)设定模拟区域的边界条件,选取粒子间作用势模型。 (3)设定系统所有粒子的初始位置和初始速度。 (4)计算粒子间的相互作用力和势能,以及各个粒子的位置和速度。 (5)待体系达到平衡,统计获得体系的宏观特性。 分子动力学模拟的主要对象就是将实际物理模型抽象后的物理系统模型。因此,物理建模也是分子动力学模拟的一个重要的环节。而对于分子动力学模拟,主要还是势函数的选取,势函数是分子动力学模拟计算的核心。这是因为分子动力学模拟主要是计算分子间作用力,计算粒子的势能、位置及速度都离不开势函数的作用。系统中粒子初始位置的设定最好与实际模拟模型相符,这样可以使系统尽快达到平衡。另外,粒子的初始速度也最好与实际系统中分子的速度相当,这样可以减少计算机的模拟时间。 要想求解粒子的运动状态就必须把运动方程离散化,离散化的方法有经典Verlet算法、蛙跳算法(Leap-frog)、速度Veriet算法、Gear预估-校正法等。这些算法有其各自的优势,选取时可按照计算要求选择最合适的算法。 统计系统各物理量时,便又涉及到系统是选取了什么系综。只有知道了模拟系统采用的系综才能釆用相对应的统计方法更加准确,有效地进行统计计算,减少信息损失。 2. 分子动力学模拟程序流程 具体到分子动力学模拟程序的具体流程,主要包括: (1)设定和模拟相关的参数。 (2)模拟体系初始化。 (3)计算粒子间的作用力。 (4)求解运动方程。 (5)循环计算,待稳定后输出结果。 分子动力学模拟程序流程图如2.3所示。
分子动力学 分子动力学方法是一种计算机模拟实验方法，是研究凝聚态系统的有力工具。该技术不仅可以得到原子的运动轨迹，还可以观察到原子运动过程中各种微观细节。它是对理论计算和实验的有力补充。 分子动力学总是假定原子的运动服从某种确定的描述，这种描叙可以牛顿方程、拉格朗日方程或哈密顿方程所确定的描述，也就是说原子的运动和确定的轨迹联系在一起。在忽略核子的量子效应和Born-Oppenheimer绝热近似下，分子动力学的这一种假设是可行的[1]。所谓绝热近似也就是要求在分子动力学过程中的每一瞬间电子都处于原子结构的基态。要进行分子动力学模拟就必须知道原子间的相互作用势。在分子动力学模拟中，我们一般采用经验势来代替原子间的相互作用势，如Lennard-Jones势、Mores势、EAM原子嵌入势、F-S多体势。然而采用经验势必然丢失了局域电子结构之间存在的强相关作用信息，即不能得到原子动力学过程中的电子性质[1]。 事实上，分子动力学就是模拟原子系统的趋衡过程。实际上，分子动力学方法就是确定某一描述与初始条件、边值关系的数值解。我们假定系统经过M步长之后达到稳定，而这
一稳定状态正是我们所求的。 １、分子动力学的算法分析 首先，我们假定我们研究的系统服从 Newton 方程所确定的描述，即： )(1 )(.. t F m t r = （1） 式中r(t)表征原子在t 时刻的位置矢量 F(t)表征原子在t 时刻所受到的力，它与所有原子的位置矢有关 m 表征原子的质量。 如果我们给定初始条件，即方程（１）的定解条件r(0)和v(0)，那么方程（１）的解就可以确定。６０年代中期发展了大量的分子动力学算法，如两步差分算法[2]、预测-校正算法 [3] 、中心差分算法[4]、蛙跳算法[5]等等。为了方便导出它们， 我们以Euler 一步法[6] 来讨论之。我们令)()(.. t r t v =（表征粒子 的速度），则有： ) ()()(1 )()(... . t v t r t F m t r t v === （２） 记??? ? ??????=? ? ? ???=)()(1)()()()(. t v t F m t f t r t v t w （３）
分子动力学模拟 分子动力学就是一门结合物理,数学与化学的综合技术。分子动力学就是一套分子模拟方法,该方法主要就是依靠牛顿力学来模拟分子体系的运动,以在由分子体系的不同状态构成的系统中抽取样本,从而计算体系的构型积分,并以构型积分的结果为基础进一步计算体系的热力学量与其她宏观性质。 这门技术的发展进程就是: 1980年:恒压条件下的动力学方法(Andersenの方法、Parrinello-Rahman法) 1983年:非平衡态动力学方法(Gillan and Dixon) 1984年:恒温条件下的动力学方法(能势‐フーバーの方法) 1985年:第一原理分子动力学法(→カー?パリネロ法) 1991年:巨正则系综的分子动力学方法(Cagin and Pettit)、 最新的巨正则系综,即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触,此时系统不仅同热源有能量交换,而且可以同粒子源有粒子的交换,最后达到平衡,这种系综称巨正则系综。 进行分子动力学模拟的第一步就是确定起始构型,一个能量较低的起始构型就是进行分子模拟的基础,一般分子的其实构型主要就是来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度,这一速度就是根据玻尔兹曼分布随机生成,由于速度的分布符合玻尔兹曼统计,因此在这个阶段,体系的温度就是恒定的。另外,在随机生成各个原子的运动速度之后须进行调整,使得体系总体在各个方向上的动量之与为零,即保证体系没有平动位移。 由上一步确定的分子组建平衡相,在构建平衡相的时候会对构型、温度等参数加以监控。 进入生产相之后体系中的分子与分子中的原子开始根据初始速度运动,可以想象其间会发生吸引、排斥乃至碰撞,这时就根据牛顿力学与预先给定的粒子间相互作用势来对各个例子的运动轨迹进行计算,在这个过程中,体系总能量不变,但分子内部势能与动能不断相互转化,从而体系的温度也不断变化,在整个过程中,体系会遍历势能面上的各个点,计算的样本正就是在这个过程中抽取的。 用抽样所得体系的各个状态计算当时体系的势能,进而计算构型积分。 作用势的选择与动力学计算的关系极为密切,选择不同的作用势,体系的势能面会有不同的形状,动力学计算所得的分子运动与分子内部运动的轨迹也会不同,进而影响到抽样的结果与抽样结果的势能计算,在计算宏观体积与微观成分关系的时候主要采用刚球模型的二体势,计算系统能量,熵等关系时早期多采用Lennard-Jones、morse势等双体势模型,对于金属计算,主要采用morse势,但就是由于通过实验拟合的对势容易导致柯西关系,与实验不符,因此在后来的模拟中有人提出采用EAM等多体势模型,或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但就是对于二体势模型,多体势往往缺乏明确的表达式,参量很多,模拟收敛速度很慢,给应用带来很大困难,因此在一般应用中,通过第一性原理计算结果拟合势函数的L-J,morse等势模型的应用仍非常广泛。 分子动力学计算的基本思想就是赋予分子体系初始运动状态之后,利用分子的自然运动在相空间中抽取样本进行统计计算,时间步长就就是抽样的间隔,因而时间步长的选取对动力学模拟非常重要。太长的时间步长会造成分子间的激烈碰撞,体系数据溢出;太短的时间步长会降低模拟过程搜索相空间的能力,因此一般选取的时间步长为体系各个自由度中最短运动周期的十分之一。但就是通常情况下,体系各自由度中运动周期最短的就是各个化学键的振动,而这种运动对计算某些宏观性质并不产生影响,因此就产生了屏蔽分子内部振动或其她无关运动的约束动力学,约束动力学可以有效地增长分子动力学模拟时间步长,提高搜索相空间的能
分子动力学方法模拟基本步骤 1.第一步 即模型的设定，也就是势函数的选取。势函数的研究和物理系统上对物质的描述研究息息相关。最早是硬球势，即小于临界值时无穷大，大于等于临界值时为零。常用的是LJ势函数，还有EAM势函数，不同的物质状态描述用不同的势函数。 模型势函数一旦确定，就可以根据物理学规律求得模拟中的守恒量。 2 第二步 给定初始条件。运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。如：verlet算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步的坐标或者一组零时刻的速度值。 一般意思上讲系统的初始条件不可能知道，实际上也不需要精确选择代求系统的初始条件，因为模拟实践足够长时，系统就会忘掉初始条件。当然，合理的初始条件可以加快系统趋于平衡的时间和步伐，获得好的精度。 常用的初始条件可以选择为：令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到；令初始位置随机的偏离差分划分网格的格子上，初始速度为零；令初始位置随机的偏离差分划分网格的格子上，初始速度也是从玻尔兹曼分布随机抽样得到。 第三步 趋于平衡计算。在边界条件和初始条件给定后就可以解运动方程，进行分子动力学模拟。但这样计算出的系统是不会具有所要求的系统的能量，并且这个状态本身也还不是一个平衡态。 为使得系统平衡，模拟中设计一个趋衡过程，即在这个过程中，我们增加或者从系统中移出能量，直到持续给出确定的能量值。我们称这时的系统已经达到平衡。这段达到平衡的时间成为驰豫时间。 分子动力学中，时间步长的大小选择十分重要，决定了模拟所需要的时间。为了减小误差，步长要小，但小了系统模拟的驰豫时间就长了。因此根据经验选择适当的步长。如，对一个具有几百个氩气Ar分子的体系，lj势函数，发现取h为0.01量级，可以得到很好的相图。这里选择的h是没有量纲的，实际上这样选择的h对应的时间在10-14s的量级呢。如果模拟1000步，系统达到平衡，驰豫时间只有10-11s。 第四步 宏观物理量的计算。实际计算宏观的物理量往往是在模拟的最后揭短进行的。它是沿相空间轨迹求平均来计算得到的（时间平均代替系综平均）
分子动力学模拟 分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。 这门技术的发展进程是： 1980年：恒压条件下的动力学方法（Andersenの方法、Parrinello-Rahman法） 1983年：非平衡态动力学方法(Gillan and Dixon) 1984年：恒温条件下的动力学方法（能势‐フーバーの方法） 1985年：第一原理分子动力学法（→カー?パリネロ法） 1991年：巨正则系综的分子动力学方法(Cagin and Pettit). 最新的巨正则系综，即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触，此时系统不仅同热源有能量交换，而且可以同粒子源有粒子的交换，最后达到平衡，这种系综称巨正则系综。 进行分子动力学模拟的第一步是确定起始构型，一个能量较低的起始构型是进行分子模拟的基础，一般分子的其实构型主要是来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度，这一速度是根据玻尔兹曼分布随机生成，由于速度的分布符合玻尔兹曼统计，因此在这个阶段，体系的温度是恒定的。另外，在随机生成各个原子的运动速度之后须进行调整，使得体系总体在各个方向上的动量之和为零，即保证体系没有平动位移。 由上一步确定的分子组建平衡相，在构建平衡相的时候会对构型、温度等参数加以监控。 进入生产相之后体系中的分子和分子中的原子开始根据初始速度运动，可以想象其间会发生吸引、排斥乃至碰撞，这时就根据牛顿力学和预先给定的粒子间相互作用势来对各个例子的运动轨迹进行计算，在这个过程中，体系总能量不变，但分子内部势能和动能不断相互转化，从而体系的温度也不断变化，在整个过程中，体系会遍历势能面上的各个点，计算的样本正是在这个过程中抽取的。 用抽样所得体系的各个状态计算当时体系的势能，进而计算构型积分。 作用势的选择与动力学计算的关系极为密切，选择不同的作用势，体系的势能面会有不同的形状，动力学计算所得的分子运动和分子内部运动的轨迹也会不同，进而影响到抽样的结果和抽样结果的势能计算，在计算宏观体积和微观成分关系的时候主要采用刚球模型的二体势，计算系统能量，熵等关系时早期多采用Lennard-Jones、morse势等双体势模型，对于金属计算，主要采用morse势，但是由于通过实验拟合的对势容易导致柯西关系，与实验不符，因此在后来的模拟中有人提出采用EAM等多体势模型，或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但是对于二体势模型，多体势往往缺乏明确的表达式，参量很多，模拟收敛速度很慢，给应用带来很大困难，因此在一般应用中，通过第一性原理计算结果拟合势函数的L-J，morse等势模型的应用仍非常广泛。 分子动力学计算的基本思想是赋予分子体系初始运动状态之后，利用分子的自然运动在相空间中抽取样本进行统计计算，时间步长就是抽样的间隔，因而时间步长的选取对动力学模拟非常重要。太长的时间步长会造成分子间的激烈碰撞，体系数据溢出；太短的时间步长会降低模拟过程搜索相空间的能力，因此一般选取的时间步长为体系各个自由度中最短运动周期的十分之一。但是通常情况下，体系各自由度中运动周期最短的是各个化学键的振动，而这种运动对计算某些宏观性质并不产生影响，因此就产生了屏蔽分子内部振动或其他无关运动的约束动力学，约束动力学可以有效地增长分子动力学模拟时间步长，提高搜索相空间的能
vasp做分子动力学的好处，由于vasp是近些年开发的比较成熟的软件，在做电子scf速度方面有较好的优势。 缺点：可选系综太少。 尽管如此，对于大多数有关分子动力学的任务还是可以胜任的。 主要使用的系综是NVT和NVE。 下面我将对主要参数进行介绍！ 一般做分子动力学的时候都需要较多原子，一般都超过100个。 当原子数多的时候，k点实际就需要较少了。有的时候用一个k点就行，不过这都需要严格的测试。通常超过200个原子的时候，用一个k点，即Gamma点就可以了。 INCAR： EDIFF 一般来说，用1E-4或者1E-5都可以，这个参数只是对第一个离子步的自洽影响大一些，对于长时间的分子动力学的模拟，精度小一点也无所谓，但不能太小。 IBRION=0 分子动力学模拟 IALGO=48 一般用48，对于原子数较多，这个优化方式较好。 NSW=1000 多少个时间步长。 POTIM=3 时间步长,单位fs,通常1到3. ISIF=2 计算外界的压力. NBLOCK= 1 多少个时间步长，写一次CONTCAR，CHG和CHGCAR，PCDAT. KBLOCK=50 NBLOCK*KBLOCK个步长写一次XDATCAR. ISMEAR=-1 费米迪拉克分布. SIGMA =0.05 单位:电子伏 NELMIN=8 一般用6到8,最小的电子scf数.太少的话,收敛的不好. LREAL=A APACO=10 径向分布函数距离,单位是埃. NPACO=200 径向分布函数插的点数. LCHARG=F 尽量不写电荷密度,否则CHG文件太大. TEBEG=300 初始温度. TEEND=300 终态温度。不设的话，等于TEBEG. SMASS -3 NVE ensemble；-1 用来做模拟退火；大于0 NVT 系综。 ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// 1)收敛判据的选择 结构弛豫的判据一般有两种选择：能量和力。这两者是相关的，理想情况下，能量收敛到基态，力也应该是收敛到平衡态的。但是数值计算过程上的差异导致以二者为判据的收敛速度差异很大，力收敛速度绝大部分情况下都慢于能量收敛速度。这是因为力的计算是在能量的基础上进行的，能量对坐标的一阶导数得到力。计算量的增大和误差的传递导致力收敛慢。 到底是以能量为收敛判据，还是以力为收敛判据呢？关心能量的人，觉得以能量
分子动力学攻略 此文为dddc_redsnow发表于biolover上的关于分子动力学的系列原创文章，相当经典与精彩，特此将系列文章整合，一起转载，望学习动力学的新手们共同学习，提高进步，在此特向dddc_redsnow本人表示感谢。 动力学系列之一（gromacs，重发） 在老何的鼓励下，发一下我的gromacs上手手册（我带人时用的，基本半天可以学会gromcas） ###################################################### # Process protein files step by step # ###################################################### pdb2gmx -f 2th_cap.pdb -o 2th_cap.gro -p 2th_cap.top -ignh -ter nedit 2th_cap.top editconf -f 2th_cap.gro -o 2th_cap_box.gro -d 1.5 genbox -cp 2th_cap_box.gro -cs -p 2th_cap.top -o 2th_cap_water.gro make_ndx -f 2th_cap_water.gro -o 2th_cap.ndx genpr -f 2th_cap_water.gro -n 2th_cap.ndx -o 2th_cap_All.itp genpr -f 2th_cap_water.gro -n 2th_cap.ndx -o 2th_cap_M.itp genpr -f 2th_cap_water.gro -n 2th_cap.ndx -o 2th_cap_C.itp nedit Flavo.itp grompp -f em.mdp -c 2th_cap_water.gro -p 2th_cap.top -o prepare.tpr genion -s prepare.tpr -o 2th_cap_water_ion.gro -np 1 -pq 1 ##################################################### # Minimize step by step # # 1. minimization fixing whole protein # # 2. minimization fixing maincharin of protein # # 3. minimization fixing Ca of protein # # 4. minimization without fix # ##################################################### grompp -np 4 -f em.mdp -c 2th_cap_water_ion.gro -p 2th_cap.top -o minimize_water.tpr mpirun -np 4 mdrun -nice 0 -s minimize_water.tpr -o minimize_water.trr -c minimize_water.gro -e minimize_water.edr -g minimize_water.log & grompp -np 4 -f em.mdp -c minimize_water.gro -p 2th_cap.top -o minimize_sidechain.tpr mpirun -np 4 mdrun -nice 0 -s minimize_sidechain.tpr -o minimize_sidechain.trr -c minimize_sidechain.gro -e minimize_sidechain.edr -g minimize_sidechain.log & grompp -np 4 -f em.mdp -c minimize_sidechain.gro -p 2th_cap.top -o minimize_sidechain_ex.tpr mpirun -np 4 mdrun -nice 0 -s minimize_sidechain_ex.tpr -o minimize_sidechain_ex.trr -c minimize_sidechain_ex.gro -e minimize_sidechain_ex.edr minimize_sidechain_ex.log & grompp -np 4 -f em.mdp -c minimize_sidechain_ex.gro -p 2th_cap.top -o minimize_all.tpr mpirun -np 4 mdrun -nice 0 -s minimize_all.tpr -o minimize_all.trr -c minimize_all.gro -e minimize_allx.edr -g minimize_all.log&
分子动力学在材料科学中的应用 摘要：本文综述了几种常见条件下的分子动力学模拟方法以及分子动力学模拟的最新发展趋势.介绍用分子动力学模拟方法研究固休的休相结构,表面问题,界面问题以及薄膜形成过程等方面的研究成果。 关键词：分子动力学; 计算机模拟; 材料科学 1引言 分子动力学（Molecular Dyanmica,简称MD）用于计算以固体、液体、气体为模型的单个分子运动，它是探索各种现象本质和某些新规律的一种强有力的计算机模拟方法，具有沟通宏观特性与微观结构的作用，对于许多在理论分析和实验观察上难以理解的现象可以做出一定的解释［1］。MD方法不要求模型过分简化，可以基于分子（原子、离子）的排列和运动的模拟结果直接计算求和以实现宏观现象中的数值估算。可以直接模拟许多宏观现象，取得和实验相符合或可以比较的结果，还可以提供微观结构、运动以及它们和体系宏观性质之间关系的极其明确的图象［2］。MD以其不带近似、跟踪粒子轨迹、模拟结果准确［3］，而倍受研究者的关注，在物理、化学、材料、摩擦学等学科及纳米机械加工中得到广泛而成功的应用。本文主要评述MD方法在材料科学中的应用. 目前在材料微观结构的研究中，由于实验条件的限制，使得许多重要的微观结构的信息难以得到，如，对于由液态金属快速凝固的非晶转变过程，其微观结构的瞬时变化根本无法用实验仪器去测量。理论分析、实验测定及模拟计算已成为现代材料科学研究的3种主要方法[2]。20世纪90年代以来，由于计算机科学和技术的飞速发展，模拟计算的地位日渐突显。计算机模拟可以提供实验上尚无法获得或很难获得的信息。虽然计算机模拟不能完全取代实验，但可以用来指导实验的进行，从而促进理论和实践的发展，所以有必要对这一领域进行介绍。 2 分子动力学基本原理 分子动力学将连续介质看成由N个原子或分子组成的粒子系统，各粒子之间的作用力可以通过量子力学势能函数求导得出，忽略量子效应后，运用经典牛顿力学建立系统粒子运动数学模型，通过数值求解得到粒子在相空间的运动轨迹，然后由统计物理学原理得出该系统相应的宏观动态、静态特性。图1所示是MD
分子动力学模拟
分子动力学模拟 分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。 这门技术的发展进程是： 1980年：恒压条件下的动力学方法（Andersenの方法、Parrinello-Rahman法）1983年：非平衡态动力学方法(Gillan and Dixon) 1984年：恒温条件下的动力学方法（能势‐フーバーの方法） 1985年：第一原理分子动力学法（→カー?パリネロ法） 1991年：巨正则系综的分子动力学方法(Cagin and Pettit). 最新的巨正则系综，即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触，此时系统不仅同热源有能量交换，而且可以同粒子源有粒子的交换，最后达到平衡，这种系综称巨正则系综。 进行分子动力学模拟的第一步是确定起始构型，一个能量较低的起始构型是进行分子模拟的基础，一般分子的其实构型主要是来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度，这一速度是根据玻尔兹曼分布随机生成，由于速度的分布符合玻尔兹曼统计，因此在这个阶段，体系的温度是恒定的。另外，在随机生成各个原子的运动速度之后须进行调整，使得体系总体在各个方向上的动量之和为零，即保证体系没有平动位移。由上一步确定的分子组建平衡相，在构建平衡相的时候会对构型、温度等参数加以监控。
进入生产相之后体系中的分子和分子中的原子开始根据初始速度运动，可以想象其间会发生吸引、排斥乃至碰撞，这时就根据牛顿力学和预先给定的粒子间相互作用势来对各个例子的运动轨迹进行计算，在这个过程中，体系总能量不变，但分子内部势能和动能不断相互转化，从而体系的温度也不断变化，在整个过程中，体系会遍历势能面上的各个点，计算的样本正是在这个过程中抽取的。 用抽样所得体系的各个状态计算当时体系的势能，进而计算构型积分。 作用势的选择与动力学计算的关系极为密切，选择不同的作用势，体系的势能面会有不同的形状，动力学计算所得的分子运动和分子内部运动的轨迹也会不同，进而影响到抽样的结果和抽样结果的势能计算，在计算宏观体积和微观成分关系的时候主要采用刚球模型的二体势，计算系统能量，熵等关系时早期多采用Lennard-Jones、morse势等双体势模型，对于金属计算，主要采用morse 势，但是由于通过实验拟合的对势容易导致柯西关系，与实验不符，因此在后来的模拟中有人提出采用EAM等多体势模型，或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但是对于二体势模型，多体势往往缺乏明确的表达式，参量很多，模拟收敛速度很慢，给应用带来很大困难，因此在一般应用中，通过第一性原理计算结果拟合势函数的L-J，morse等势模型的应用仍非常广泛。 分子动力学计算的基本思想是赋予分子体系初始运动状态之后，利用分子的自然运动在相空间中抽取样本进行统计计算，时间步长就是抽样的间隔，因而时间步长的选取对动力学模拟非常重要。太长的时间步长会造成分子间的激烈碰撞，体系数据溢出；太短的时间步长会降低模拟过程搜索相空间的能力，因此
《材料计算设计基础》 学号： 流水号： 姓名： 完成日期：
分子动力学模拟及其在材料中的研究进展 摘要：本文综述了分子动力学模拟技术的发展，介绍了分子动力学的分类、运动方程的求解、初始条件和边界条件的选取、平衡系综及其控制、感兴趣量的提取以及分子动力学模拟在材料中的研究进展。 关键词：分子动力学模拟平衡态系综金属材料感兴趣量径向分布函数 引言 科学工作者在长期的科学研究实践中发现，当实验研究方法不能满足研究工作的需求时，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；尽管计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从而促进理论和实验的发展。特别是在材料形成过程中许多与原子有关的微观细节，在实验中基本上是无法获得的，而在计算机模拟中即可以方便地得到。这种优点使分子动力学模拟在金属材料研究中显得非常有吸引力。 分子动力学MD (Molecular Dynamics)模拟就是用计算机方法来表示统计力学，作为实验的一个辅助手段。MD模拟就是对于原子核和电子所构成的多体系统，求解运动方程(如牛顿方程、哈密顿方程或拉格朗日方程)，其中每一个原子核被视为在全部其它原子核和电子作用下运动，通过分析系统中各粒子的受力情况，用经典或量子的方法求解系统中各粒子在某时刻的位置和速度，以确定粒子的运动状态，进而计算系统的结构和性质。该模拟技术主要涉及粒子运动的动力学问题,与蒙特卡罗模拟方法(简称MC)相比,分子动力学是一种“确定性方法”, 它所计算的是时间平均,而MC进行的是系综平均。然而按照统计力学各态历经假设,时间平均等价于系综平均。因此,两种方法严格的比较计算能给出几乎相同的结果。 经典的分子动力学方法是Alder等于1957年提出并首先在“硬球”液体模型下应用，发现了由Kirkwood在1939年根据统计力学预言的“刚性球组成的集合系统会发生有液相到结晶相的转变”。后来人们称这种相变为Alder相变。Rahman
分子动力学模拟基础知识 ? Molecular Dynamics Simulation o MD: Theoretical Background Newtonian Mechanics and Numerical Integration The Liouville Operator Formalism to Generating MD Integration Schemes o Case Study 1: An MD Code for the Lennard-Jones Fluid Introduction The Code, mdlj.c o Case Study 2: Static Properties of the Lennard-Jones Fluid (Case Study 4 in F&S) o Case Study 3: Dynamical Properties: The Self-Diffusion Coefficient ? Ensembles o Molecular Dynamics at Constant Temperature Velocity Scaling: Isokinetics and the Berendsen Thermostat Stochastic NVT Thermostats: Andersen, Langevin, and Dissipative Particle Dynamics The Nosé-Hoover Chain Molecular Dynamics at Constant Pressure: The Berendsen Barostat Molecular Dynamics Simulation We saw that the Metropolis Monte Carlo simulation technique generates a sequence of states with appropriate probabilities for computing ensemble averages (Eq. 1). Generating states probabilitistically is not the only way to explore phase space. The idea behind the Molecular Dynamics (MD) technique is that we can observe our dynamical system explore phase space by solving all particle equations of motion . We treat the particles as classical objects that, at least at this stage of the course, obey Newtonian mechanics. Not only does this in principle provide us with a properly weighted sequence of states over which we can compute ensemble averages, it additionally gives us time-resolved information, something that Metropolis Monte Carlo cannot provide. The ``ensemble averages'' computed in traditional MD simulations are in practice time averages : (99) The ergodic hypothesis partially requires that the measurement time, , i , in the system. The price we pay for this extra information is that we must at least access if not store particle velocities in addition to positions, and we must compute interparticle forces in addition to potential energy. We will introduce and explore MD in this section.
分子动力学(MD) 1 分子动力学（MD）基础 1.1 MD分类 1.2 MD简介 1.3 MD适用范围 2 分子动力学运动方程数值求解 2.1 基础知识 2.1.1 运动方程 2.1.2 空间描述 2.1.3 最小作用量原理 2.1.4 拉格朗日(Lagrange)方程 2.1.5 哈密顿(Hamilton)方程 2.2 粒子运动方程的数值解法 2.2.1 Verlet算法 2.2.2 欧拉(Euler)预测—矫正公式 2.2.3 Gear预测—矫正方法 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 分子力场函数的构成
4.2.2 常用力场函数和分类 5 分子动力学模拟的基本步骤 5.1 设定模拟所采用的模型 5.2 给定初始条件 5.3 趋于平衡计算 5.4 宏观物理量的计算 6 平衡态分子动力学模拟 6.1 系综 6.2 微正则系综的分子动力学模拟6.3 正则系综的分子动力学模拟
1 分子动力学（MD）基础 1.1MD分类 微正则系综（VNE） 正则系综（VNP） 平衡态MD 等温等压系综（NPT） 经典MD 等焓等压系综（NPH） 巨正则系综（ＶＴμ） 非平衡态ＭＤ 量子MD 1.2分子动力学(MD)简介 分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法。分子动力学方法为确定性模拟方法，广泛地用于研究经典的多粒子体系的研究中,是按该体系内部的内禀动力学规律来计算并确定位形的转变。 分子动力学方法是通过建立一组分子的运动方程，并通过直接对系统中的一个个分子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计计算方法得到多体系统的静态和动态特性, 从而得到系统的宏观性质。 在分子动力学中，粒子的运动行为是通过经典的Newton运动方程所描述。系统的所有粒子服从经典力学的运动规律，它的动力学方程就是从经典力学的运动方程——拉格朗日(lagrange)方程和哈密顿(Hamilton)方程导出。 1.3适用范围 原则上，分子动力学方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分子，也可以是其它的微观粒子。 实际上，分子动力学模拟方法和随机模拟方法一样都面临着两个基本限制：
分子动力学模拟 The Standardization Office was revised on the afternoon of December 13, 2020
分子动力学模拟 分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。 这门技术的发展进程是： 1980年：恒压条件下的动力学方法（Andersenの方法、Parrinello-Rahman法）1983年：非平衡态动力学方法(Gillan and Dixon) 1984年：恒温条件下的动力学方法（能势‐フーバーの方法） 1985年：第一原理分子动力学法（→カー?パリネロ法） 1991年：巨正则系综的分子动力学方法(Cagin and Pettit). 最新的巨正则系综，即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触，此时系统不仅同热源有能量交换，而且可以同粒子源有粒子的交换，最后达到平衡，这种系综称巨正则系综。 进行分子动力学模拟的第一步是确定起始构型，一个能量较低的起始构型是进行分子模拟的基础，一般分子的其实构型主要是来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度，这一速度是根据玻尔兹曼分布随机生成，由于速度的分布符合玻尔兹曼统计，因此在这个阶段，体系的温度是恒定的。另外，在随机生成各个原子的运动速度之后须进行调整，使得体系总体在各个方向上的动量之和为零，即保证体系没有平动位移。由上一步确定的分子组建平衡相，在构建平衡相的时候会对构型、温度等参数加以监控。
1、分子动力学简介： 分子动力学方法是一种计算机模拟的实验方法，是研究凝聚态系统的有力工具。该技术不仅可以得到原子的运动轨迹，还可以观察到原子运动过程中各种微观细节。它是对理论计算和实验的有力补充。广泛应用于材料科学、生物物理和药物设计等。经典MD模拟，其系统规模在一般的计算机上也可达到数万个原子，模拟时间为纳秒量级。 分子动力学总是假定原子的运动服从某种确定的描述，这种描叙可以牛顿方程、拉格朗日方程或哈密顿方程所确定的描述，也就是说原子的运动和确定的轨迹联系在一起。在忽略核子的量子效应和绝热近似（Born-Oppenheimer）下，分子动力学的这一种假设是可行的。所谓绝热近似也就是要求在分子动力学过程中的每一瞬间电子都处于原子结构的基态。要进行分子动力学模拟就必须知道原子间的相互作用势。 在分子动力学模拟中，我们一般采用经验势来代替原子间的相互作用势，如Lennard-Jones势、Mores势、EAM原子嵌入势、F-S多体势。然而采用经验势必然丢失了局域电子结构之间存在的强相互作用的信息，即不能得到原子动力学过程中的电子性质。
2、分子模拟的三步法和大致分类 三步法： 第一步：建模。包括几何建模，物理建模，化学建模，力学建模。初始条件的设定，这里要从微观和宏观两个方面进行考虑。 第二步：过程。这里就是体现所谓分子动力学特点的地方。包括对运动方程的积分的有效算法。对实际的过程的模拟算法。关键是分清楚平衡和非平衡，静态和动态以及准静态情况。 第三步：分析。这里是做学问的关键。你需要从以上的计算的结果中提取年需要的特征，说明你的问题的实质和结果。因此关键是统计、平均、定义、计算。比如温度、体积、压力、应力等宏观量和微观过程量是怎么联系的。 大致分类： 2.1电子模拟（量化计算，DFT） 量子化学计算 一般处理几个到几十个原子 常见软件：GAUSSIAN，NWCHEM等 密度泛函（DFT） 可以算到上百个原子 常见软件：V ASP 2.2分子模拟（分子动力学，蒙特卡洛） 2.2.1分子级别的模拟 以分子的运动为主要模拟对象。采用经验性的分子间作用函数