一．分子动力学简介
随着纳米科技的到来，许多新的学科产生了，例如纳米电子学、纳米生物学、纳米材料学、纳米机械学等。人们的注意力逐渐从宏观物体转向小尺度及相应的器件，其中微机械系统(mieromachine)或称微型机电系统(mieroe一eetro一meeh耐ealsystem，MEMs)尤其取得了成功，并正被拓展应用于各种工业过程。由图1可知，分子动力学正是处于nm尺度下的研究方法。
图1.不同模拟方法所对应的空间和时间尺度
1957年Alder和Wainwright[1]开创了分子动力学(Moleeularnynamies，MD)方法，之后经过多位科学家的努力，拓展了分子动力学方法的理论、技术及应用领域，尤其是在20世纪80年代由Andersen等[2]先后完成的恒温、定压分子动力学方法，标志着分子动力学方法的科研应用进入了一个新阶段。
分子动力学方法是研究纳米尺度物理现象的重要手段。随着越来越多的材料原子间作用势函数被精确描述并经过实验验证、计算机硬件水平的快速更新以及高效率新算法的提出，分子动力学模拟被广泛应用于纳米尺度力学行为和纳米材料力学性能的研究。
在纳米尺度下，材料由离散的原子排列而成，由于比表面积大、表面效应明显，材料的力学性能和力学行为将与宏观材料迥异。基于连续性假设的宏观连续介质理论在研究材料的损伤演化、失效过程时，往往在时间和空间上将原子尺度的缺陷进行平均化处理，但这种处理仅适用于大量缺陷分布在材料中计算区域的情形，而对许多细微观材料和力学实验观测到的现象都无法解释，如疲劳与蠕变过程中的位错模式、塑性变形的不均匀性、脆性断裂的统计本质、尺寸效应等。因此，连续介质理论显然难以准确求解纳米尺度的力学问题。同时，如果直接从第一原理出发进行计算，除了类氢原子以外其他材料的薛定愕方程求解难度都太大，而且局域密度泛函近似理论并不是总能满足实际问题的需要。另一方面，材料本身在空间、时间和能量等方面存在藕合和脱祸现象[3,4]，直接从头开始的量子力学计算难以很好地应用到几百个原子以下的计算规模中，无法达到一般纳米材料和器件的模拟要求。此外，由于实验条件控制的困难和合成、制备方式不同，各种纳米材料力学性能的有关实验结果分散性较大甚至相反，以至于目前难以通过纳米力学实验得到普遍适用的定量力学规律。
鉴于理论和实验上的困难，通过分子动力学方法模拟纳米尺度的力学性能和行为来探索纳米尺度的一般规律，是进行纳米力学研究的有效方法。分子动力学最早用于热动力学和物理化学，计算不同物理系统如固体、气体、液体的整体或平均热化学性能。1957年Alder 首次提出并采用分子动力学方法分析刚性球系统的固液相变问题取得成功，此后，分子动力
学开始逐渐应用于材料领域。随着上世纪80年代计算机硬件水平的提高和各种描述原子间作用的势函数的提出，分子动力学模拟日益活跃。通过分子动力学模拟不仅可以深入了解复杂的机制，发现本质上崭新的现象，而且可以定量模拟真实固体中所发生的过程，是诸如表面结构和扩散中的动力学和稳定性等许多结果的唯一来源[5-9]。在EAM理论逐渐成熟和Baskes实验室[10-13]、Ackland实验室[14-18]等精确测出大量常用材料的EAM参数以后，分子动力学方法在模拟材料的物理性能和现象方面逐渐显示出强大的计算能力和较高的精度，大量的模拟尤其是固体结构、位错运动、表面界面现象、力学性能、变形机制和流体中的电泳、电渗透流等都得到了理想的结果。只要能将基于物理的模型建立起来，通过分子动力学计算就可以揭示出物理现象的本质，逐渐被广泛应用于凝聚态物理、材料学、力学、生物学、微电子学和微纳米加工等领域。
晶体是由大量的原子有序排列而成，材料的强度来源于原子间的相互作用，塑性来源于原子间的相互运动。因此，直接从原子尺度上对材料的微观力学行为进行研究显得非常有必要。分子动力学模拟技术既能得到原子的运动轨迹，还能像做实验一样观察。对于平衡系统，可以在一个分子动力学观察时间(ObservationTime)内做时间平均来计算一个物理量的统计平均值，对于一个非平衡系统过程，只要发生在一个分子动力学观察时间内(一般为1一10ps)的物理现象也可以用分子动力学计算进行直接模拟。特别是许多与原子有关的微观细节，在实际实验中无法获得，而在计算机模拟中可以方便的得到。这种优点使得分子动力学方法广泛运用于材料科学与工程中，如材料设计、断裂分析等。
近年来，利用计算机模拟技术研究材料的力学性能日益成为人们感兴趣的课题。由于计算机处理速度的迅速提高，计算机模拟已经和实验观察、理论分析并列成为本世纪科学研究的三种方法[19]。计算机模拟的数据(从模型中得来)可以用来比较、验证各种近似理论;同时，计算机模拟方法还可用来对实验和模型进行比较，从而提供了评估一个模型正确与否的手段。计算机模拟方法还有一个优点就是可以沟通理论和实验。某些量或行为可能是无法或难以在实验中测量的，而用计算机模拟方法，这些量可以被精确的计算出来[20]。分子动力学模拟方法更以其建模简单、模拟结果准确的特征而倍受研究者的关注。
分子动力学模拟(MolecularDynamicsSimulation)用于计算以固体、液体、气体为模型的单个分子的运动情况，是一种联系微观世界与宏观世界的强有力的计算机模拟方法侧。
二．分子动力学基本原理
分子动力学方法是纳米计算力学的主要手段。它对多经典粒子系统的运动方程组进行数值积分，得到相空间轨道，进而研究系统的平衡热力学性质、结构动力学性质和非平衡传输性质等。在分子动力学模拟中，将原子视为质点或忽略内部结构的固体球。首先建立粒子系统的几何模型，通过描述原子间相互作用的势函数，根据给定的边界条件、初始条件求出系统中每一时刻单个粒子或原子的能量和所受到的力，代入牛顿动力学方程组求解原子的位置和速度，得到系统在相空间的运动轨迹。对足够大数量的粒子在足够长时间的结果进行统平均，则可以得到类似于宏观意义上的物理量和力学量。
早期研究孤立系统的保守系综分子动力学基于两个基本假设[21]:
a.粒子是相互作用的质点，运动由位置矢量和速度矢量来描述。粒子间的相互作用取决于彼此的空间位置。
b.系统无质量交换。即模拟过程中系统的原子数不变。
保守系综分子动力学常用于模拟能量守恒的孤立绝热系统。但是要准确模拟实际纳米材料和器件的表面、界面等现象，真实反映粒子系统受外界约束情况下的物理行为，以及进行跨尺度的祸合模拟，保守系综分子动力学模拟明显是远远不够的。在此基础上，研究者提出了耗散系统与周围介质进行能量交换的不同理论和算法，对粒子系统进行温度、压力、粒子数、体积等不同控制，适合于不同系综的分子动力学模拟。
分子动力学模拟中假设系统所有粒子的运动遵循经典牛顿运动定律，且忽略原子背景电子云的量子效应，原子间的相互作用满足叠加原理。可见，分子动力学方法是一种对广义牛顿运动方程的近似的数值积分方法。
1.下面分别描述分子动力学方法的基本方程。
(1)拉格朗日（Lagrangian）运动方程
对于N自由度的由互相作用的质点构成的运动系统，拉格朗日方程为：
式中q为广义坐标，指定质点的空间位置。q 为质点位置对时间的导数。在笛
卡儿坐标系下，由N个原子构成的模拟系统的拉格朗日方程可写为：
式中，为原子i的位置矢量，在三维空间n=N/3。对于互不影响的粒子系统，可取拉格朗日量：
式中m i，为粒子i的质量，L即为系统的动能，即系统所有原子的动能之和。如果考虑到原子间的互相作用，L可改一记为
式右端的两项分别表示系统的动能和势能。代入拉格朗日方程可得系统牛顿运动方程:
式中F i即为原了i所受的内力，即由于系统中其他原子的作用而在原子i上体现出的合力。由牛顿运动方程建立线性微分方程组，给定初始条件(初始位置、初始速度)，求解该封
闭方程，可得到确定解，求出原子运动的轨迹，即任意时刻原子的位置r i(t)和速度.
r
i (t)。
分子动力学模拟中的经典拉格朗日方程常用来计算原子的运动，得到单个原子的运动轨迹，描述原子系统的运动过程，以及反映在原子系统整体特征下的位移、变形、缺陷等。
(2)哈密顿(Hamiltonian)运动方程：
如果采用广义坐标系和动量的形式来描述粒子系统的运动，则可以得到哈密顿形式的运动方程，求出多粒子系统的状态和演化过程。
在笛卡儿坐标系中，对含n个粒子的保守系统，取系统哈密顿量为：
式中P i
，为粒子i 的动量，P i =m
i r
i ，哈密顿正则方程为:
如果给定系统的初始状态(初始位置和速度)，则可以求解线性微分方程组，得到系统原子的运动轨迹，即任意时刻的原子位置和动量，由统计平均得到系统的热力学表征。
在分子动力学模拟的基本方程中，拉格朗日运动方程更适合于求解原子系统运动的过程(求解原子速度和位置)，并能施加外部荷载如外力、约束、边界条件等，而保守统的哈密顿正则方程更适合于求解系统的动力学演化过程和热力学状态，如温度、热流动。在实际的分子动力学模拟中，拉格朗日方程通常与外部约束和边界条件一起，构成特定原子系统模型的初值问题;而Hamiltonian 运动方程往往在保守系统的基础上，通过改变系统状态变量，与外部环境进行能量交换，构成不同系综的分子动力学模拟。
2.分子动力学原子间作用势函数
分子动力学方法是通过原子间的相互作用势，按照经典牛顿运动定律求出原子运动轨迹及其演化过程。分子动力学计算的关键是原子间势函数的选取，它决定着计算的工作量以及计算模型与真实系统的近似程度，直接影响到模拟结果的成功与否。由于物质系统的复杂性以及原子间相互作用类型的不同，很难得到满足各种不同体系和物质的一般性而又精度较高的势函数。所以针对不同的物质体系人们陆续发展了大量的经验和半经验的势函数[22]。
从分子动力学模拟的基本方程可以看出，分子动力学模拟即求解偏微分方程，
求解的精度关键在于势函数U(r 1,r 2,…r n )描述原子间相互作用的准确性。
多原子系统的势函数可以表示成: 式中U m 为m 体势，即原子势能受m 阶多体效应的影响，其中U 1为系统原子受外力场(如重力)的影响项。在通常的计算中，由于系统内原子受原子间相互作用的影响远大于单个原子本身受外力场的影响，同时为了减少计算量，一般忽略外力场的影响U 1和三阶以上的多体效应U m (m>3)，这就是应用于分子动力学模拟的对势模型 (PairPotential)。对不同的固体材料，将其原子三阶以上多体效应的影响作为修正项记入对势模型的二阶作用项中，就形成了各种不同的多体势函数。 （1）对偶式 对偶势理论认为任何两个原子结合键的强度不因为它周围其它原子的存在而受到影响。 1.Lennard 一Jones(L 一J)势函数 在对势模型的典型代表、应用广泛的Lennard 一Jones 模型[23-24]中，取原子间作用的势函数方程为:
式中σ—原子与零点势的距离;
ε—最小势能处的能量;
m,n—调整系数，一般(m,n)取值为(12, 6), (10, 5), C8, 4)
方程右端的第一项描述原子间的排斥作用，第二项描述原子间的吸引作用。|r y|=|r i-r j|，即原子间的距离。当原了间距离r ij=σ时，原子势能Ф(r ij)=0;r ij>σ时，Ф(r ij)<0;r ij<σ时，Ф(r ij)>0。ε表征原子间的吸引/排斥强度，对原子中的一个从平衡位置r0移到无穷远处所作的功，对于不同物质的原子,ε各不相同，通常为10-19J或ev量级。能量最低时刻的r值，即r0，表示原子间的平衡位置，此时原子间作用力为零。对上述方程求极值，可得
62σ≈1.22σ
r0=
将上式代入式，可得Lconard一JoneS 模型中原了间作用力为：
2.Morse势函数
1929年，Morse注意到双原子分子的振动谱的量子力学问题可用指数形式的势函数解析地解决，并发现计算结果与实验一致[22]。于是他提出如下形式的势函数:
式中D—结合能系数;
α—势能曲线梯度系数;
r0—分子间作用力为零时的原子间距。
原子间作用力分别描述为:
同样，等式右边第一项描述原子核间的短程排斥作用，第二项描述原子核对电
子云的吸引作用。
3.Born一Mayer势函数
Bom和Mayer[25]估计碱金属离子之间的排斥项可用指数形式表示，于是提出如下形式的势函数:
Bom一Mayer排斥项等价于Morse势函数的第一项，第三项表示偶极子一偶极子
相互作用，第四项表示偶极子一四极子相互作用。A ij反映了泡利不相容效应。该势
函数考虑了离子之间的库仑力，所以它能很好地应用于碱金属卤化物以及碱土金属
卤化物等各种离子型化合物。
4.Jonson势函数
利用对偶势人们能很好地拟合FCC金属的弹性常数但用于BCC金属却误差较大。后来，考虑了原子体积效应，Johnson在对势项的基础上提出了用于BCC金属的Johnson势函数[26]:
.(2)多体势
20世纪80年代陆续发展出许多考虑多体相互作用的新的势函数(关于中心对称的)。主要有Daw和Bakes在1983年提出的嵌入原子势，Norskov和Lang在1980年发展的有效介质理论，凝胶模型，Finnis一Sinclair多体势等。其中前三种势函数的理论基础是一致的，为我们在这里主要介绍嵌入原子势(EAM)和F一S多体势。
在模拟固体和复杂的分子结构时，要考虑原子间的多体效应，即势函数方程中的高阶项。但多体势项中Φ≥3n的引入将导致参数过多，一方面这些参数的求解和实验拟合比较困难，另一方面过多的参数将带来计算上的困难。从上世纪八十年代开始，科学家开始研究在对势的基础上进行多体效应修正，即考虑三阶以上的多体作用，将其作为一个修正项引入到对势模型中，使势函数能更精确描述固体的原子间作用。根据这一思路，研究者先后提出了一系列新的多体势函数，如能较精确地描述C，Si原子间作用的Tersoff势[27]，描述C和烃类化合物的Brenne:势[28]、适于描述体心立方晶格的Finnis一Sinclair[29]:势等等。这些势函数的共同点是在对势模型的原子间吸引项上增加一个修正参数，即将势函数写为:
式中ΦR(r ij)和ΦA(r ij)分别是对势模型中原子距离r ij=|r ij|的原子i、j之间的排斥作用项和吸引作用项。εij表征其他的邻近空间原子对当前原子对的多体效应。
1.镶嵌原子势(EAM)
镶嵌原子势(EAM)[10-13]由Daw和BaskeS于1984年率先提出。它在对势模型的基础上，假设将每个原子镶嵌在由所有其他近邻原子产生的背景电了云中，由电子云密度和将原子“镶嵌”到该电子云中所需的能量来衡量原子所处背景电子云的特征和背景电子石对势能的影响。由于在电子云密度的计算中采用球面平均，镶嵌原子法特别适用于电子云基本呈球对称分布的金属原子。这种方法单独求解电子百对原子核的多体效应，而不是在原子间对势吸引项上采取参数修正的拟合方法，结果更为精确。
依据上述原理，在镶嵌原子法中，将系统势能写为：
式中右边第一项即为传统的对势项，，第二项为多体修正项。Φ（r ij）表示原子核i、j之间的排斥作用：
z i （r ij ）和z j （r ij ）为两原子核中的有效电荷数。由于受到屏蔽作用，有效电荷数随着原子间距离r ij =|r ij |的增大而衰减：
式中Z i 为原子的价电子数，α为待定常数，由材料的弹性模量确定。
F(ρi )是把第i 个原子核“镶嵌”到密度为ρi 的背景电子云中所需的镶嵌能，可以写为： F(ρi )=D ρi In ρi
D 是由材料决定的参数，ρi 即为除第i 个原子以外其他所有原子在原子i 处所产生的电子云密度的线性叠加：
ρi=()r ij
i j j ∑Φ≠
Фj （r ij ）是离原子i 处距离为r ij 的j 原子在原子i 处产生的电子云密度。
在EAM 模型提出并得到逐步完善后，特别是计算结果得到大量实验验证后，原子间作用的精确描述问题逐渐转移到测试和求解EAM 势参数上来。
针对原始的镶嵌原子求解电子云密度时采用球面平均而不适用于极性原子键结构，并对部分材料预测不准，暴露出对堆垛层错计算过大的缺点，Baskes 研究组对EAM 进行了部分修正，并于上世纪90年代先后测出了元素周期表上常用元素的EAM 参数，使得EAM 模型得到更广泛和方便的应用。 Doyama 等从第一原理出发给出了求解立方晶格金属原子镶嵌势参数的方法，并得到实验验证。在他们的求解中，将势函数方程中的对势项取为：
电子云密度假设为：
以上两式中r c1和r c2分别为计算原子间作用项和电子云密度项的截断半径。这样，立方晶格金属晶体势函数的求解就集中在求解由材料物理属性如结合能、玻恩稳定性、弹性常数、空位形成能和堆垛能等决定的5个参数A 1、A 2、c 1、c 2、D 上。
2.Finnis 一Sinelair 势[25]
F 一S 势和EAM 势的理论基础不同，它是建立在紧束缚(TB)近似的基础上的。F 一S 势函数可以写成如下形式:
V 和Φ均为经验拟合的对势函数，式(2一23)中第一项表示核一核之间的相互作用，第二项不再是嵌入能，而表示多体相互作用的关联能，Ф可以理解为二分量紧束缚近似下的积分平方求和。Finnis 一Sinclair 模型最初用来描述具有BCC 结构的纯金属的原子间相互作用，Ackland 等人将它推广到二元合金体系。对于二元合金体系，可以用下面6个对势函数表示 V AA 、V BB 、V AB 、ΦAA 、ΦBB 、ΦAB 。Ackland 等认为V AA 、V BB 、ΦAA 、ΦBB 和纯金属的参数相同；函数ΦAB 取ΦAA 和ΦBB 的几何平均，这和紧束缚近似下的积分求和的解释是一致的。因此在已知纯金属势参数的情况下，仅有对势项屹。需要通过拟合而得到。
Ackland给出上述模型的函数:
H(x)是阶跃函数。与EAM相比，F一S多体势形式较为简单，势参数易于拟合。此为计算量比较少，程序容易收敛，所以应用广泛。
3.分子动力学模拟模型
在给定原子间作用势函数的条件下，可以建立模拟系统的几何模型，根据模拟的实际物理情况施加初始条件和边界条件，选取时间步长，然后按照一定的动力学积分算法，通过求解牛顿运动方程组或系统Hamilton正则运动方程组，计算系统中原子在每一步的状态，包括每个原子的位置、速度和原子势，并由原子势计算每个原子上的作用力，根据原子力给出下一时间步的原子位置和速度。在模拟过程中对原子系统实施相应的边界条件调节，实现系统和周围介质的能量交换。对足够长时间的计算结果和足够大系统的原子计算结果进行统计平均，可以得到类似于宏观概念的系统物理量。
(1)几何模型的建立
几何模型的建立就是按照要模拟的真实晶体中原子排列给出原子系统构型。根据晶体结构分析，绝大多数金属晶体均为面心立方晶格((Face Centered Crystal ,FCC，如Au, Ag, Cu, Ni, Al,γ-Fe等)、体心立方晶格(Body Centered Crystal,BCC，如Cr, Mo, W, Nb, V,α-Fe等)和密排六方晶格(Hexagonal Closed Packed,HCP，如Mg, Cd, Zn, Be等)这三种典型的紧密结构，单个晶胞的原子排列结构如下图。空间方向上依次排列多个单胞，规定相应的品格尺寸，就构成了要模拟的单晶纳米材料，同时可以根据模拟需要在模型中引入晶格缺陷，如表面、裂纹等。如果要构造多晶材料，则需要按照一定的分布形式如Voronoi几何构造方法构筑多个单晶晶粒，形成晶界。
图.常见金属晶体结构
(2)分子动力学模拟中的有限差分算法
分子动力学模拟的主要目的是计算任意时刻相互作用的原子系统中各原子的运动轨迹，即求解动力学系统初值问题，理论上可以采用有限差分法求解系统的牛顿运动方程。分子动力学模拟中，粒子的运动方程是一组常微分方程组，求解的基本原理都是利用有限差分法，以一定的时间步长对方程沿时间轴逐步积分，但分子动力学对计算时间有一些特殊要求。首先由于计算分子间作用力所占的工作量极大，所以凡是要求在一个时间步中重复计算分子间作用力的算法是不适宜的。第二方面的要求涉及到解的精度和稳定性，这要求算法的阶数至少应大于3。适用于分子动力学模拟的算法有Verlet算法、Leap-frog算法、Gear算法、Euler算法等
1) Verlet算法
Verlet方法基于有限差分思想，由初始时刻的状态值及其时间导数得到下一时
刻的状态值，并由此逐步计算出整个物理过程。在分子动力学计算中，对原子的状
态进行Taylor级数展开:
两式相加可得:
由中两式相减可得原子速度:
根据系统的牛顿运动方程可得原子加速度:
这就是Verlet算法的基本形式。
2) Gear算法
为了求解非线性常微分方程，Gear提出了基于预测一校正积分方法的Gear算法[30]。Gear0预测一校正法通过Taylor级数展开，由当前时刻原子的位置、速度、加速度以及其他高阶导数计算下一时刻原子的位置、速度、加速度，即预报值，然后根据原子位置的预报值计算原子间作用力和实际加速度值，通过比较实际加速度值和预报加速度值，进行误差校正，得到下一时刻的原子位置、速度、加速度以及其他高阶导数。将t+?t时刻的位置、速度和加速度对时间t作泰勒展开可得预测步为：
由新给出的位置r p(t十?t)可计算t+?t时刻的力F(t + ?t)和加速度αr(t+?t),根据加速度αr 评价预测的加速度αp (t + ?t)的误差为:
将此差值与预测项相加得修正项为:
式中C0、C1、C2和C3为修正项系数。
在分子动力学模拟中，除上述三种常用的算法外，还有很多其它一些用于特殊要求情况下的算法。Swope提出的Velocity-Verlet算法可以同时给出位置、速度与加速度，并且不牺牲精度。这种算法的优点是给出了显式速度项，计算量适中。Beeman提出的Beeman算法运用了更精确的速度表达式，所以能更精确地计算系统动能。然而，由于它的表达式很复杂，所以计算量很大[31]。Satoh就各种分子动力学算法的稳定性所作的对比研究表明，最优越的算法是Velocity-Verlet算法，之后是Leap-frog算法，Beeman算法等。Gear算法比这些算法都要差些，因其更复杂，且稳定性和能量波动性较差，因而只是常见于早期的分子动力学模拟以及简单系统的计算。
参考文献：
[ 1 ] Alder B J, Wainwright T E. Phase Transition for Hard Sphere System. The Journal of Chemical P坤sics. 1957,27(5):1208-1209.
[2] Stillinger F H, Weber T A. Phys Rev B. 1985,31:52-62.
[3]冯瑞，金国均.凝聚态物理中的基本概念.物理学进展.2000, 20 (1) :1-21.
[4」胡壮麟，王鲁红.电子和原子层次材料行为的计算机模拟.材料研究学报.1998,
12 (1):1一19.
[5〕梁海戈.纳米铜力学行为的分子动力学模拟〔博士论文].合肥:中国科技大学.2001.
[6〕吴恒安，倪向贵，王宇等.金属纳米棒弯曲力学行为的分子动力学模拟.物理学报.2002, 51(7):1412-1415.
[7〕吴恒安，王秀喜，梁海戈等.金属纳米丝力学行为研究进展.金属学报.2002,
38 (9):903-907.
[8〕吴恒安，王秀喜，倪向贵等.金属纳米杆弹性模量的直接原子计算.金属学
报.2002, 38 (11):1219-1222.
[9〕吴恒安.纳米尺度下结构和材料力学行为的分子动力学模拟研究[博士论文〕.
合肥:中国科技大学.2002.
[10] Daw BS, Baskes MI. Embedded-atom method: Derivation and application to
impurities surfaces and other defects in metals. P勿sical Review B.1984, 29:
6443-6453.
[川Baskes MI. Application of the semiempirical potential for silicon. method to covalent materials: A Physical Review Letters.1987, 59: 2666-2669.
[12] Baskes MI, Nelson JS, Wright AF. Semiempirical embedded-atom potential for
silicon and germanium. Physical Review B.1989, 40: 6085-6100.
[13] Baskes MI. Modified embedded-atom potentials for cubic materials and impurities.
Physical Review B.1992, 46: 2727-2742.
[ 14] Ackland GJ, Thetford R. An improved N-body semi-empirical model for body-centered cubic transition metals. Philosophical MagazineA.1987, 56(1):
15-30.
[15] Ackland GJ, Ticy GSVitek V,et al.Simple N-body potentials for the noble metals
and nickel. Philosophical Magazine A.1987, 56(3): 735-764.
[16] Ackland GT. Theoretical study of titaniwn surfaces and defects with a new
many-body potential. Philosophical Magazine A.1992, 66(6): 917-932.
[17] Ackland GJ, Mendelev MI, Srolovitz DJ, et al. Development of phosphorus ties in
an interatomic potential for alpha-iron. Journal of P勿sics: Condensed Matter.2004, 16: 2629-2642.
[18] Ackland GJ, Bacon DJ, Calder AF, et al. Computer simulation of point defect
properties in dilute Fe-Cu alloy using a many-body interatomic potential.
Philosophical Magazine A.1997, 75(3): 713-732.
[19]罗旋.β-SiC表面及A1/SiC界面结构的分子动力学模拟〔博士论文〕.黑龙江:
哈尔滨工业大学.1997.
[20] Cotterill R M, Doyama M. Energy and atomic configuration of complete and
dissociated dislocations[J]. PhysRev.1966,145:465.
[21] Allen M D, Tildesley D J, Computer Simulation of Liquid, Oxford: Clarendon Press.
1987: 256-317.
[22]陈强，草红红，黄海波.分子动力学中势函数研究.天津理工学院学
报.2004, 20 (2):101一105.
[23]Jones JE. On the determination of molecular fields. I .From the variation of the
viscosity of a gas with tempe Proceedings of Royal Society A. 1924, 106: 441-462.
[24]Jones JE. On the determination of molecular fields. II .From the equation of state of
a gas. Proceedings of Royal Society A. 1924, 106: 463-477.
[25]黄海波.L1o-TiAl中角度相关势和Ni3A1中点缺陷的分子动力学研究[硕士论
文〕.北京:北京航空航天大学.2003.
[26]Wilson W D, Johnson R A. Non-Central-Force Model of LiH: Phonon Dispersion
Curves and He Migration. Phys. Rev. B. 1970,8:3510-3514.
[27]Tersoff J. Empirical interatomic potential for carbon with applications to
amorphous carbon. P勿sical Review Letters.l988, 61(25): 2879-2882.
[28]Rosenblum I, Adler J, Brandon S. Multi-processor molecular dynamics using the
Brenner potential: parallelization of an implicit mufti-body optional. International Journal of Modem P场sics C.1999, 10(1): 189-203.
[29] Finnis MW, Sinclair JE. A simple empirical N-body potential for transition metals.
Philosophical Magazine A. 1984, 50(1): 545-554.
[30]樊康旗.基于两态稳定性的超高密度数据存储〔硕士论文〕.陕西:西安电子
科技大学.2005.
[31]文玉华，朱如曾，周富信等.分子动力学模拟的主要技术.力学进
展.2003, 33 (1):65-73.
《装备制造技术》 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所示。
物理计算与设计报告书 院（系）名称： 学生姓名： 专业名称： 班级： 时间： 金属铝分子动力学模拟
摘要：分子动力学模拟，是指对于原子核和电子所构成的多体系统，用计算机模拟 原子核的运动过程，并从而计算系统的结构和性质，其中每一原子核被视为在全部其它 原子核和电子所提供的经验势场作用下按牛顿定律运动。我们用c语言编写程序，VMD 动画演示得到原子在拉伸过程中的变化。在控制温度不变的情况下，得到了金属铝分子 的动力学模拟过程。通过不断拉伸，趋衡铝分子，计算其势能，力,速度，观察每次拉 伸过程中以及拉伸后铝原子的排列，得到金属铝的运动细节，从而更加利于我们了解铝 的性质。 结论：原子两端的拉力与原子势能的变化曲线基本一致。原子间断层以滑层方式断 裂。 关键词：铝分子，分子动力学，c语言，势能 1 引言 人们很早就知道材料的力学性能随尺度发生变化尺度减小, 材料中缺陷存在的几率降低, 材料的强度提高同时尺度的变化可能导致材料内在变形竞争机制的改变, 例如多晶材料晶粒粒径在微米级以上时, 强度主要受位错强化机制控制, 而粒径进入纳米级后, 材料的变形主要来源于晶界滑移等机制原子尺度下, 微观效应占主导地位, 材料的理化、力学性能表现出与宏观不同、甚至相反的特性。Brenner发现金属单晶晶须拉伸强度与晶须直径呈反比,Fleck在微米级细铜丝的扭转试验中观察到尺寸效应纳米电机系统(NEMS)的出现同迫切要求了解纳米尺度下材料的力学行为, 当前从实验上较难获得详细的信息, 而分子动力学模拟可以提供相关细节. 分子动力学通过直接模拟原子的运动过程, 使我们能够详细了解模拟对象的演化发展历史分子动力学模拟的一个关键在于原子势函数的选取原子势早期一般采用简单的对势, 但对势无法正确描述弹性常数, 其结果不理想世纪年代提出的镶嵌原子法、有效介质理论更客观地反映了原子间多体作用的本质, 可得到较合理的结果.认为体系总能量为
分子动力学方法模拟基本步骤 1.第一步 即模型的设定，也就是势函数的选取。势函数的研究和物理系统上对物质的描述研究息息相关。最早是硬球势，即小于临界值时无穷大，大于等于临界值时为零。常用的是LJ势函数，还有EAM势函数，不同的物质状态描述用不同的势函数。 模型势函数一旦确定，就可以根据物理学规律求得模拟中的守恒量。 2 第二步 给定初始条件。运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。如：verlet算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步的坐标或者一组零时刻的速度值。 一般意思上讲系统的初始条件不可能知道，实际上也不需要精确选择代求系统的初始条件，因为模拟实践足够长时，系统就会忘掉初始条件。当然，合理的初始条件可以加快系统趋于平衡的时间和步伐，获得好的精度。 常用的初始条件可以选择为：令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到；令初始位置随机的偏离差分划分网格的格子上，初始速度为零；令初始位置随机的偏离差分划分网格的格子上，初始速度也是从玻尔兹曼分布随机抽样得到。 第三步 趋于平衡计算。在边界条件和初始条件给定后就可以解运动方程，进行分子动力学模拟。但这样计算出的系统是不会具有所要求的系统的能量，并且这个状态本身也还不是一个平衡态。 为使得系统平衡，模拟中设计一个趋衡过程，即在这个过程中，我们增加或者从系统中移出能量，直到持续给出确定的能量值。我们称这时的系统已经达到平衡。这段达到平衡的时间成为驰豫时间。 分子动力学中，时间步长的大小选择十分重要，决定了模拟所需要的时间。为了减小误差，步长要小，但小了系统模拟的驰豫时间就长了。因此根据经验选择适当的步长。如，对一个具有几百个氩气Ar分子的体系，lj势函数，发现取h为0.01量级，可以得到很好的相图。这里选择的h是没有量纲的，实际上这样选择的h对应的时间在10-14s的量级呢。如果模拟1000步，系统达到平衡，驰豫时间只有10-11s。 第四步 宏观物理量的计算。实际计算宏观的物理量往往是在模拟的最后揭短进行的。它是沿相空间轨迹求平均来计算得到的（时间平均代替系综平均）
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)收敛判据的选择 结构弛豫的判据一般有两种选择：能量和力。这两者是相关的，理想情况下，能量收敛到基态，力也应该是收敛到平衡态的。但是数值计算过程上的差异导致以二者为判据的收敛速度差异很大，力收敛速度绝大部分情况下都慢于能量收敛速度。这是因为力的计算是在能量的基础上进行的，能量对坐标的一阶导数得到力。计算量的增大和误差的传递导致力收敛慢。 到底是以能量为收敛判据，还是以力为收敛判据呢？关心能量的人，觉得以能量
分子动力学模拟 分子动力学就是一门结合物理,数学与化学的综合技术。分子动力学就是一套分子模拟方法,该方法主要就是依靠牛顿力学来模拟分子体系的运动,以在由分子体系的不同状态构成的系统中抽取样本,从而计算体系的构型积分,并以构型积分的结果为基础进一步计算体系的热力学量与其她宏观性质。 这门技术的发展进程就是: 1980年:恒压条件下的动力学方法(Andersenの方法、Parrinello-Rahman法) 1983年:非平衡态动力学方法(Gillan and Dixon) 1984年:恒温条件下的动力学方法(能势‐フーバーの方法) 1985年:第一原理分子动力学法(→カー?パリネロ法) 1991年:巨正则系综的分子动力学方法(Cagin and Pettit)、 最新的巨正则系综,即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触,此时系统不仅同热源有能量交换,而且可以同粒子源有粒子的交换,最后达到平衡,这种系综称巨正则系综。 进行分子动力学模拟的第一步就是确定起始构型,一个能量较低的起始构型就是进行分子模拟的基础,一般分子的其实构型主要就是来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度,这一速度就是根据玻尔兹曼分布随机生成,由于速度的分布符合玻尔兹曼统计,因此在这个阶段,体系的温度就是恒定的。另外,在随机生成各个原子的运动速度之后须进行调整,使得体系总体在各个方向上的动量之与为零,即保证体系没有平动位移。 由上一步确定的分子组建平衡相,在构建平衡相的时候会对构型、温度等参数加以监控。 进入生产相之后体系中的分子与分子中的原子开始根据初始速度运动,可以想象其间会发生吸引、排斥乃至碰撞,这时就根据牛顿力学与预先给定的粒子间相互作用势来对各个例子的运动轨迹进行计算,在这个过程中,体系总能量不变,但分子内部势能与动能不断相互转化,从而体系的温度也不断变化,在整个过程中,体系会遍历势能面上的各个点,计算的样本正就是在这个过程中抽取的。 用抽样所得体系的各个状态计算当时体系的势能,进而计算构型积分。 作用势的选择与动力学计算的关系极为密切,选择不同的作用势,体系的势能面会有不同的形状,动力学计算所得的分子运动与分子内部运动的轨迹也会不同,进而影响到抽样的结果与抽样结果的势能计算,在计算宏观体积与微观成分关系的时候主要采用刚球模型的二体势,计算系统能量,熵等关系时早期多采用Lennard-Jones、morse势等双体势模型,对于金属计算,主要采用morse势,但就是由于通过实验拟合的对势容易导致柯西关系,与实验不符,因此在后来的模拟中有人提出采用EAM等多体势模型,或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但就是对于二体势模型,多体势往往缺乏明确的表达式,参量很多,模拟收敛速度很慢,给应用带来很大困难,因此在一般应用中,通过第一性原理计算结果拟合势函数的L-J,morse等势模型的应用仍非常广泛。 分子动力学计算的基本思想就是赋予分子体系初始运动状态之后,利用分子的自然运动在相空间中抽取样本进行统计计算,时间步长就就是抽样的间隔,因而时间步长的选取对动力学模拟非常重要。太长的时间步长会造成分子间的激烈碰撞,体系数据溢出;太短的时间步长会降低模拟过程搜索相空间的能力,因此一般选取的时间步长为体系各个自由度中最短运动周期的十分之一。但就是通常情况下,体系各自由度中运动周期最短的就是各个化学键的振动,而这种运动对计算某些宏观性质并不产生影响,因此就产生了屏蔽分子内部振动或其她无关运动的约束动力学,约束动力学可以有效地增长分子动力学模拟时间步长,提高搜索相空间的能
分子动力学模拟 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引言 分子动力学（Molecular Dyanmica,简称MD）用于计算以固体、液体、气体为模型的单个分子运动，它是探索各种现象本质和某些新规律的一种强有力的计算机模拟方法，具有沟通宏观特性与微观结构的作用，对于许多在理论分析和实验观察上难以理解的现象可以做出一定的解释［1］。MD方法不要求模型过分简化，可以基于分子（原子、离子）的排列和运动的模拟结果直接计算求和以实现宏观现象中的数值估算。可以直接模拟许多宏观现象，取得和实验相符合或可以比较的结果，还可以提供微观结构、运动以及它们和体系宏观性质之间关系的极其明确的图象［2］。MD以其不带近似、跟踪粒子轨迹、模拟结果准确［3］，而倍受研究者的关注，在物理、化学、材料、摩擦学等学科及纳米机械加工中得到广泛而成功的应用。本文主要评述MD方法在材料科学中的应用. 目前在材料微观结构的研究中，由于实验条件的限制，使得许多重要的微观结构的信息难以得到，如，对于由液态金属快速凝固的非晶转变过程，其微观结构的瞬时变化根本无法用实验仪器去测量。理论分析、实验测定及模拟计算已成为现代材料科学研究的3种主要方法[2]。20世纪90年代以来，由于计算机科学和技术的飞速发展，模拟计算的地位日渐突显。计算机模拟可以提供实验上尚无法获得或很难获得的信息。虽然计算机模拟不能完全取代实验，但可以用来指导实验的进行，从而促进理论和实践的发展，所以有必要对这一领域进行介绍。 2 分子动力学基本原理 分子动力学将连续介质看成由N个原子或分子组成的粒子系统，各粒子之间的作用力可以通过量子力学势能函数求导得出，忽略量子效应后，运用经典牛顿力学建立系统粒子运动数学模型，通过数值求解得到粒子在相空间的运动轨迹，然后由统计物理学原理得出该系统相应的宏观动态、静态特性。图1所示是MD
第六章 分子动力学模拟 Molecular Dynamics –MD 6.1引言 分子动力学模拟方法是在牛顿力学的理论框架下，根据体系内分子之间的相互作用势，获得每个原子随时间运动的轨迹，通过系综平均，可以得到感兴趣的与结构和动力学性质有关的物理量，如：平均原子坐标，平均能量、平均温度及原子运动的自相关函数等。这些物理量是通过对每个原子的运动轨迹，即微观量求平均而得到的宏观量，因此可以与实验观测量进行比较。 用计算机模拟方法在向空间采样方法有两种： （1） 随机采样 MC （2） 确定性方法MD 以上讲过的MC （Monte Carlo ）采样方法就是随机方法，与随机方法不同，确定性方法是按照动力学规律使系统在相空间运动。分子动力学模型就是一种确定性方法。它的基本出发点是从一个完全确定的物理模型出发，通过解牛顿运动方程而得到原子运动的轨迹。我们感兴趣的可测量的客观物理量可以通过相空间的采样求系综平均而得到。在多态历经假设成立的情况下，系综平均与长时间平均是相同的。 ? ∞ →∞==τ τ τ0 1 ))(),((lim dt t p t q A A A 系综 其中q,p 为t 的函数。A 表示系综平均，∞A 表示无穷长时间平均。因模拟时间总是有限的。对耦分子体系，当模拟时间大于分子的弛豫时间时，有限观测时间可以变成为无穷长的。 当弛豫模拟〉τt ，模拟t 可认为∞，因物理上的∞是不可能的。 6.2基本原理 1.动力学方程 基本动力学方程包括在经典力学（CM ）框架下的牛顿方程和在量子动力学（QM ）框架下的薛定谔方程。在常温下，经典的牛顿方程对研究生物分子体系的结构和动力学性质已经足够了，因为这时体系的量子效应并不十分重要。但是，对研究包含隧道效应的反应时间问题时，量子效应十分明显，这时就必须用QM 方程来模拟体系的量子动力学性质。
分子动力学模拟 分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。 这门技术的发展进程是： 1980年：恒压条件下的动力学方法（Andersenの方法、Parrinello-Rahman法） 1983年：非平衡态动力学方法(Gillan and Dixon) 1984年：恒温条件下的动力学方法（能势‐フーバーの方法） 1985年：第一原理分子动力学法（→カー?パリネロ法） 1991年：巨正则系综的分子动力学方法(Cagin and Pettit). 最新的巨正则系综，即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触，此时系统不仅同热源有能量交换，而且可以同粒子源有粒子的交换，最后达到平衡，这种系综称巨正则系综。 进行分子动力学模拟的第一步是确定起始构型，一个能量较低的起始构型是进行分子模拟的基础，一般分子的其实构型主要是来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度，这一速度是根据玻尔兹曼分布随机生成，由于速度的分布符合玻尔兹曼统计，因此在这个阶段，体系的温度是恒定的。另外，在随机生成各个原子的运动速度之后须进行调整，使得体系总体在各个方向上的动量之和为零，即保证体系没有平动位移。 由上一步确定的分子组建平衡相，在构建平衡相的时候会对构型、温度等参数加以监控。进入生产相之后体系中的分子和分子中的原子开始根据初始速度运动，可以想象其间会发生吸引、排斥乃至碰撞，这时就根据牛顿力学和预先给定的粒子间相互作用势来对各个例子的运动轨迹进行计算，在这个过程中，体系总能量不变，但分子内部势能和动能不断相互转化，从而体系的温度也不断变化，在整个过程中，体系会遍历势能面上的各个点，计算的样本正是在这个过程中抽取的。 用抽样所得体系的各个状态计算当时体系的势能，进而计算构型积分。 作用势的选择与动力学计算的关系极为密切，选择不同的作用势，体系的势能面会有不同的形状，动力学计算所得的分子运动和分子内部运动的轨迹也会不同，进而影响到抽样的结果和抽样结果的势能计算，在计算宏观体积和微观成分关系的时候主要采用刚球模型的二体势，计算系统能量，熵等关系时早期多采用Lennard-Jones、morse势等双体势模型，对于金属计算，主要采用morse势，但是由于通过实验拟合的对势容易导致柯西关系，与实验不符，因此在后来的模拟中有人提出采用EAM等多体势模型，或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但是对于二体势模型，多体势往往缺乏明确的表达式，参量很多，模拟收敛速度很慢，给应用带来很大困难，因此在一般应用中，通过第一性原理计算结果拟合势函数的L-J，morse等势模型的应用仍非常广泛。 分子动力学计算的基本思想是赋予分子体系初始运动状态之后，利用分子的自然运动在相空间中抽取样本进行统计计算，时间步长就是抽样的间隔，因而时间步长的选取对动力学模拟非常重要。太长的时间步长会造成分子间的激烈碰撞，体系数据溢出；太短的时间步长会降低模拟过程搜索相空间的能力，因此一般选取的时间步长为体系各个自由度中最短运动周期的十分之一。但是通常情况下，体系各自由度中运动周期最短的是各个化学键的振动，而这种运动对计算某些宏观性质并不产生影响，因此就产生了屏蔽分子内部振动或其他无关运动的约束动力学，约束动力学可以有效地增长分子动力学模拟时间步长，提高搜索相空间的能力。
分子动力学方法 一、引言 计算机模拟中的另一类确定性模拟方法，即统计物理中的所谓合于动力学方法（Molecular Dynamics Method）。这种方法是按该体系内部的内禀动力学规律来计算并确定位形的转变。它首先需要建立一组分子的运动方程，并通过直接对系统中的一个个分子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计计算方法得到多体系统的静态和动态特性，从而得到系统的宏观性质。在这样的处理过程中我们可以看出：MD方法中不存在任何随机因素。在MD方法处理过程中方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学。每个分子运动的内禀动力学是用理论力学上的哈密顿量或者拉格朗日量来描述，也可以直接用牛顿运动方程来描述。确定性方法是实现Boltzman的统计力学途径。这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。但是使用该方法的程序较复杂，讨算量大，占内存也多、本节将介绍分子动力学方法及其应用。 原则上，MD方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分子，也可以是其他的微观粒子。 实际上，MD模拟方法和随机模拟方法一样都面临着两个基本限制：一个是有限观测时间的限制；另一个是有限系统大小的限制。通常人们感兴趣的是体系在热力学极限下（即粒子数日趋于无穷时）的性质。但是计算机模拟允许的体系大小要比热力学极限小得多，因此可能会出现有限尺寸效应。为了减小有限尺寸效应，人们往往引入周期性、全反射、漫反射等边界条件。当然边界条件的引入显然会影响体系的某些性质。 对于MD方法，向然的系综是微正则系综，这时能量是运动常量。然而，当我们想要研究温度和（或）压力是运动常量的系统时，系统不再是封闭的。例如当温度为常量的系统可以认为系统是放置在一个热俗中。当然，在MD方法中我们只是在想像中将系统放入热浴中。实际上，在模拟计算中具体所采取的做法是对一些自由度加以约束。例如在恒温体系的情况下，体系的平均动能是一个不变量。这时我们可以设计一个算法，使平均动能被约束在一个给定值上。由于这个约束，我们并不是在真正处理一个正则系综，而实际上仅仅是复制了这个系综的位形部分。只要这一约束不破坏从一个状态到另一个状态的马尔科夫特性，这种做法就是正确的。不过其动力学性质可能会受到这一约束的影响。 自五十年代中期开始，MD方法得到了广泛的应用。它与蒙特卡洛方法一起已经成为计算机模拟的重要方法。应用MD方法取得了许多重要成果，例如气体或液体的状态方程、相变问题、吸附问题等，以及非平衡过程的研究。其应用已从化学反应、生物学的蛋白质到重离子碰撞等广泛的学科研究领域。 二、分子运动方程的数值求解 采用MD方法时，必须对一组分于运动微分方程做数值求解。从计算数学的角度来看，这个求解是一个初值问题。实际上计算数学为了求解这种问题己经发展了许多的算法，但并不是所有的这些算法都可以用来解决物理问题。下面我们先以一个一维谐振子为例，来看一下如何用计算机数值计算方法求解初值问题。一维谐振子的经典哈密顿量为 （2.1） 这里的哈密顿量（即能量）为守恒量。假定初始条件为x(p)、p(0),则它的哈密顿方程是对时间的一阶微分方程 （2.2） 现在我们要用数值积分方法计算在相空间中的运动轨迹(X(t)、p(t)) 。我们采用有限差分法，将微分方程变为有限差分方程，以便在计算机上做数值求解，并得到空间坐标和动量随时间的演化关系。首先，
分子动力学模拟方法的基本原理与应用 摘要: 介绍了分子动力学模拟的基本原理及常用的原子间相互作用势, 如Lennard-Jones势; 论述了几种常用的有限差分算法, 如Verlet算法; 说明了分子动力学模拟的几种系综及感兴趣的宏观统计量的提取。 关键词: 分子动力学模拟; 原子间相互作用势; 有限差分算法; 分子学是一门结合物理，和化学的综合技术。分子学是一套方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的量和其他宏观性质。 从统计物理学中衍生出来的分子动力学模拟方法(Molecular Dynamics Simulation, MDS) , 实践证明是一种描述纳米科技研究对象的有效方法, 得到越来越广泛的重视。所谓分子动力学模拟, 是指对于原子核和电子所构成的多体系统, 用计算机模拟原子核的运动过程, 从而计算系统的结构和性质, 其中每一个原子核被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动。它被认为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段, 称之为“计算机实验”手段, 在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。 科学工作者在长期的科学研究实践中发现，当实验研究方法不能满足研究工作的需求时，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；尽管计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从而促进理论和实验的发展。特别是在材料形成过程中许多与原子有关的微观细节，在实验中基本上是无法获得的，而在计算机模拟中即可以方便地得到。这种优点使分子动力学模拟在材料研究中显得非常有吸引力。 分子动力学模拟就是用计算机方法来表示统计力学，作为实验的一个辅助手段。分子模拟就是对于原子核和电子所构成的多体系统，求解运动方程(如牛顿方程、哈
分子动力学(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适用范围 原则上，分子动力学方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分子，也可以是其它的微观粒子。 实际上，分子动力学模拟方法和随机模拟方法一样都面临着两个基本限制：
分子动力学模拟课程小结 一．分子动力学的基本原理 在分子动力学模拟中，体系原子的一系列位移是通过对牛顿运动方程积分得到的，结果是一条运动轨迹，它表明了系统内原子的位置与速度如何随时间而发生变化。通过解牛顿第二定律的微分方程，可以获得原子的运动轨迹。方程如下： 这个方程描述了质量为m i的原子i在力Fi的作用下，位置矢量为r i时的运动方程。其中，Fi可以由势函数U的梯度给出： 系统的温度则与系统中全部原子的总动能K通过下式相联系： N是原子数，Nc是限制条件，k B是波尔兹曼常数。 二. MD模拟的积分算法 为了得到原子的运动轨迹，可以采用有限差分法来求解运动方程。有限差分法的基本思想就是将积分分成很多小步，每一小步的时间固定为δt。用有限差分解运动方程有许多方法，所有的算法都假定位置与动态性质（速度、加速度等）可以用Taylor级数展开来近似： 在分子动力学模拟中，常用的有以下的几中算法： 1.Verlet算法 运用t时刻的位置和速度及t－δt时刻的位置，计算出t＋δt时刻的位置： 两式相加并忽略高阶项，可以得到： 速度可以通过以下方法得到： 用t＋δt时刻与t－δt时刻的位置差除以2δt：
同理，半时间步t+δt时刻的速度也可以算： Verlet算法执行简单明了，存储要求适度，但缺点是位置r(t＋δt)要通过小项与非常大的两项2r(t)与r(t－δt)的差相加得到，容易造成精度损失。另外，其方程式中没有显示速度项，在没有得到下一步的位置前速度项难以得到。它不是一个自启动算法：新位置必须由t时刻与前一时刻t－δt的位置得到。在t＝0时刻，只有一组位置，所以必须通过其它方法得到t－δt的位置。一般用Taylor级数： 2.Velocity－Verlet算法 3.Leap-frog算法 为了执行Leap-frog算法，必须首先由t-0.5δt时刻的速度与t时刻的加速度计算出速度v(t+δt)，然后由方程 计算出位置r(t+δt)。T时刻的速度可以由： 得到。速度蛙跳过此t时刻的位置而得到t＋0.5δt时刻的速度值，而位置跳过速度值给出了t＋δt时刻的位置值，为计算t＋0.5δt时刻的速的作准备，依此类推。其缺点是位置与速度不同步。这意味着在位置一定时，不可能同时计算动能对总能量的贡献。 三. 分子动力学计算的时间间隔 时间间隔δt在积分算法中是一个非常重要的参数。为了充分利用CPU时间，尽量选择比较大的时间间隔，但是如果时间间隔太大，就会造成积分过程的不稳定性和不精确性。时间间隔的设置同时依赖于算法和模型的情况。模型本身给时间间隔带来的最大的限制就是最高频率的运动。由于Verlet算法要求在每个时间间隔内模型的速度和加速度保持一边，时间间隔就应该低于振动周期的八分之一到十分之一。对大多数的有机模型来讲，最高的振动频率是C-H键的伸缩振动，其振动周期的数量级为10－14s。这样，时间间隔就应该是0.5-1fs左右。如果采用受约束的SHAKE或者RATTLE算法，可以使用更长的时间间隔。如果研究对象是液态或者固态简单模型，对体系内作用模式不感兴趣，也可以采用一些更长的时间间隔，比如20fs。对离子态的材料模型，5fs左右是合适的。时间间隔必须跟选择的算法相匹配。比如，ABM4算法的时间间隔应该是Verlet算法的一半左Runge-Kutta-4
分子动力学模拟与药物研发 著：雅格布·D·德兰特，J·安德鲁·麦坎蒙，译：张浩晨 摘要：本文讨论了例如蛋白质等大分子结构的受体及其相应的小分子配体的多原子计算机模拟在药物研发中的作用。包括对于隐秘或存在变构现象的结合位点的识别，对于传统虚拟筛选方法的改进增强，以及对于小分子结合所需能量的直接预测。目前所用模拟方法的局限性，即较高的计算成本和所要求分子模拟力场的近似性，也在文中提及。随着计算机计算能力的提升和分子动力学算法的改进，计算机辅助药物设计的前景充满希望，其中分子动力学模拟技术将充当越来越重要的作用。 关键字：分子动力学模拟，计算机辅助药物设计，隐秘结合位点，变构结合位点，虚拟筛选， 自由能预测 介绍 1965年诺贝尔物理学奖获得者理查德·费曼，曾经有句名言说：“如果我们要做出一种假设来引导我们去探索生命的奥秘，那就是一切都是由原子构成的，而生命体所做的一切行为都可以理解为原子间有规律的或温和或剧烈的跃动。”过去的50年里，很多生物物理学一直致力于更好地了解原子的这种摇动和摆动的性质。微观世界的量子力学运动规律令那些只熟悉宏观动力学的人们感到惊讶。微观运动没有特定的规律，而是存在概率性；化学键的形成不是机械连接，而是通过改变以波和粒子状态存在的电子云。正如费曼所说的那样，大自然是那样的令人摸不着头脑。 理解这些荒谬而复杂但很神奇的分子运动规律无疑是与药物研发关系密切且很有帮助的。最初研究受体和配体结合的“锁和钥匙”理论认为容纳小分子配体的受体的结合域构象是如冰刻般一动不动，不存在任何的结构重排现象。现在这个观点已经基本被否定。取而代之的是一种新的观点，认为分子间的结合不仅引起构象的变化，同时也存在相结合分子间随机的摇摆运动。 软体动物的乙酰胆碱结合蛋白（AChBP）作为人类烟碱乙酰胆碱受体（AChR）的配体结合域结构和功能的替代，是无数可以用来说明原子运动重要性的例子之一（图1）。在结合于AChBP的小分子AChR激动剂的晶体结构显示，在结合处一个关键的环状结构部分闭锁（图 1 a,c）。相反，当体积较大的乙酰胆碱受体拮抗剂，如蛇的α -神经毒素结合到AchBP时的晶体结构表明，该环状结构扩大了至少10?，使得活性位点更加凸显（图1 b,c）。在研究未结合配体的样品的构象状态时显示，同时存在开放和闭合构象。伯恩等人提出未结合配体的AChBP或AChR是高度动态的蛋白质，推测蛋白构象在有选择地结合激动剂或拮抗剂后才趋于稳定。这些结合口袋构象中的任何一个都可能具有药效作用，因此也可能在药理学上有重要意义。如果这个配体结合的一般模型是正确的，则对于基于结构的药物设计的影响将是深远的。因为相似的结合原则可能也适用于许多其它体系。 分子动力学模拟 虽然这些晶体结构的研究有力地证明了蛋白质受体结合配体时其结构发生重排的特性起到重要的作用，但是得到这些结果所需要的费用和大量的劳动力，致使很多人转向去研究可以预测蛋白质运动的计算技术。不幸的是，模拟复杂的
第6章分子动力学方法 经典分子动力学方法无疑是材料，尤其是大分子体系和大体系模拟有效的方法之一。分子动力学可以用于NPT，NVE，NVT等不同系综的计算，是一种基于牛顿力学确定论的热力学计算方法。与蒙特卡罗法相比在宏观性质计算上具有更高的准确度和有效性，可以广泛应用于物理，化学，生物，材料，医学等各个领域。本章在介绍分子动力学的基本概念的基础上，简单介绍了分子动力学的基本思想，势函数分类和基本方程。然后介绍了分子动力学的常用系综和典型的NPT，NVE，NVT系综基本方程。结合材料建模中的基本简化方法和技巧，阐述了边界条件和时间积分的数值处理技巧。最后，利用统计力学的基本概念给出分子动力学的计算信息的解析方式。并且结合Materials Explore软件计算分析了CNT的几何结构稳定性。 6.1引言 分子动力学方法(Molecular Dynamics, MD)方法是一种按该体系部的禀动力学规律来计算并确定位形的变化的确定性模拟方法。首先需要在给定的外界条件下建立一组粒子的运动方程，然后通过直接对系统中的一个个粒子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计力学方法得到多体系统的静态和动态特性，从而获得系统的宏观性质。可以看出，分子动力学方法中不存在任何随机因素，这个也是分子动力学方法和后文要提到的蒙特卡洛方法的区别之一。在分子动力学方法的处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学定律(或者是拉格朗日方程)。每个分子运动的禀动力学是用理论力学上的哈密顿量或者拉格朗日函数来描述，也可以直接用牛顿运动方程来描述。确定性方法是实现玻尔兹曼的统计力学途径。这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。但是分子动力学方法的计算机程序相对蒙特卡罗较复杂，其计算成本较高。 分子动力学方法发展历史改革经历了近60年，分子动力学方法是20世纪50年代后期由Alder B J 和Wainwright T E创造发展的。Alder和Wainwright在1957年利用分子动力学模拟，发现了“刚性球组成的集合系统会发生由其液相到结晶相的相转变”，后来人们称这种相变为Alder相变。其结果表明，不具有引力的粒子系统也具有凝聚态。到20世纪70年代，产生了刚性体系的动力学方法，成功地被应用于水和氮等分子性溶液体系的处理；1972年，Less A W和 Edwards S F 等发展了该方法并扩展到了存在速度梯度，即处于非平衡状态的系统。之后，此方法被 Ginan M J等推广到了具有温度梯度的非平衡系统，从而构造并形成了所谓的非平衡分子动力学方法体系。进人20世纪80年代之后，出现了在分子部对一部分自由度施加约束条件的新的分子动力学方法，从而使分子动力学方法可适用于类似蛋白质等生物大分子的解析与设计。分子动力学方法真正作为材料科学领域的一个重要研究方法，开始于由Andersen, Parrinello和Rahman等创立恒压分子动力学方法和Nǒse 等完成恒温分子动力学方法的建立及在应用方面的成功。后来，针对势函数模型化比较困难的半导体和金属等，1985 年人们又提出了将基于密度泛函理论的电子论和分子动力学方法有机统一起来的所谓Car-Parrinello 方法，亦即第一性原理分子动力学方法。这样，分子动力学的方法进一步得到发展和完善，它不仅可以处理半导体和金属的间题，同时还可应用
第六章 分子动力学方法 6.1引言 对于一个多粒子体系的实验观测物理量的数值可以由总的平均得到。但是由于实验体系又非常大，我们不可能计算求得所有涉及到的态的物理量数值的总平均。按照产生位形变化的方法，我们有两类方法对有限的一系列态的物理量做统计平均： 第一类是随机模拟方法。它是实现Gibbs的统计力学途径。在此方法中，体系位形的转变是通过马尔科夫(Markov)过程，由随机性的演化引起的。这里的马尔科夫过程相当于是内禀动力学在概率方面的对应物。该方法可以被用到没有任何内禀动力学模型体系的模拟上。随机模拟方法计算的程序简单，占内存少，但是该方法难于处理非平衡态的问题。
另一类为确定性模拟方法，即统计物理中的所谓分子动力学方法（Molecular Dynamics Method）。这种方法广泛地用于研究经典的多粒子体系的研究中。该方法是按该体系内部的内禀动力学规律来计算并确定位形的转变。它首先需要建立一组分子的运动方程，并通过直接对系统中的一个个分子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计计算方法得到多体系统的静态和动态特性, 从而得到系统的宏观性质。因此，分子动力学模拟方法可以看作是体系在一段时间内的发展过程的模拟。在这样的处理过程中我们可以看出：分子动力学方法中不存在任何随机因素。 系统的动力学机制决定运动方程的形式： 在分子动力学方法处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学。每个分子运动的内禀动力学是用理论力学上的哈密顿量或者拉格朗日量来描述，也可以直接用牛顿运动方程来描述。这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。但是使用该方法的程序较复杂，计算量大，占内存也多。