分子动力学模拟方法的基本原理与应用
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
MD模拟结果准确与否的关键在于对系统内的原子之间相互作用势函数的选取， 总的来说，原子（或分子）之间的相互作用势的研究进展一直很缓慢，在一定程度 上制约了MD方法在实际研究中的应用，原子间的势函数的发展经历了从对势到多体 势的过程，对势认为原子之间的相互作用是两两之间的作用，与其它原子的位置无 关，而实际上，在多原子体系中，一个原子的位置不同，将影响其它原子间的有效 相互作用。所以，多体势能更准确地表示多原子体系势函数。
温度的变化为ΔT=(λ2-1)T(t)，其中，
，T 为所控制体系的温度。
（2）Berendsen 恒温槽方法：这种控温方法假设所模拟的体系与一个恒温槽连
载一起，则两者之间就可以通过热交换而使模拟体系达到恒温的目的，方法如下，
定义一个参数 λ：
，式中τ表征系Leabharlann 与恒温槽之间的热交换速率，Δt 为 MD 的时间步长，那么通过 vnew=λvold 校正即可保持体系的温度在 T0 附近震动，而参数 （通常取为~）则可用于控制这个震动幅度。
平衡态分子动力学模拟是在一定的系综瞎进行的，经常用的平衡系综是 NVT 或 NPT 系综。在这两种系综中，牵涉到控制温度和压力的几种技术，分别如下： 、控温方法：
（1）速度标定：系统的温度与动能存在如下关系：
，
式中，N 是原子数，Nc 是约束数，KB 是 Boltzmann 常数，vi 是原子 i 的速度。由于系 统的温度和动能存在这样的关系，所以一种最简单和最直观的方法是直接对速度进 行标定。这种方法的基本思想是：如果 t 时刻的温度是 T(t)，速度乘以因子λ后，
温方法类似，如下：
，式中 是一个可调参数，
在处理粒子之间的相互作用时，通常采用“最小影像”约定。这个约定是在由 无穷重复的 MD 基本模拟计算元胞中，一个粒子只与它所在的基本元胞内的另外 N-1 个（设在此元胞内有 N 个粒子） 中的每个粒子或其最邻近影像粒子发生相互作用。 实际上，这个约定就是通过满足不等式条件 rc<L/2 来截断位势（rc 为截断半径，L 是元胞的边长）通常 L 的数值应当选得很大，以避免有限尺寸效应，但这样会增大 计算量，同城采用对相互作用势的修正来近似处理。 4、势函数：
科学工作者在长期的科学研究实践中发现，当实验研究方法不能满足研究工作 的需求时，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；尽 管计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从 而促进理论和实验的发展。特别是在材料形成过程中许多与原子有关的微观细节， 在实验中基本上是无法获得的，而在计算机模拟中即可以方便地得到。这种优点使 分子动力学模拟在材料研究中显得非常有吸引力。
用 MD 方法研究压力的诱导相变和结构重构，在等压下模拟比在等体积下更容易 实现。在等压模拟下，可以通过改变模拟元胞的三个方向或一个方向的尺寸来实现 体积的变化，类似于温度控制的方法，也有许多方法用于压力控制，在实际应用中， 有以下两种较为常用的压力控制方法。
（1） Berendsen 方法：这种控压方法的基本思路与基本控制方程均与他的控
分子动力学模拟的出发点是假定粒子的运动可以用经典动力学来处理，对一个 由Ｎ个粒子构成的孤立体系，粒子的运动由牛顿运动方程决定，也就是： mid2ri/dt2=-▽iV(r1,r2,rN),式中，mi，ri 分别为第 i 个原子的质量和位置。▽ i=-?/?ri，V(r1,r2,rN)为体系所处的势。
Hockey 提出的 Leap-frog 算法是 Verlet 算法的变化，这种方法设计半时间间 隔的速度，即：r(t+δt)=r(t)+ δtv(t+δt/2)，v(t+δt/2)=v(t-δt/2)+ δtα(t)。
t 时刻的速度由下式给出：v(t)=[v(t+δt/2)+v(t-δt/2)]/2。这种算法与 Verlet 算法相比有两个优点：（1）包括显速度项；（2）收敛速度快，计算量小。这种算法 明显的缺陷是位置和速度不同步。
除了上述提及的几种方法外, 还有 Beeman 算法、Rahman 等。 3、周期性边界条件和长程力：
即使是使用现代的巨型计算机，MD 方法还是只能用于粒子数大约是几百到几千 的系统。这就引起一个问题：用这样少量的粒子，如何来模拟宏观体系为了解决这 个问题，引入了周期性边界条件[3]。采用这种方法，模拟体系实际上是由基本单元 （也称为模拟计算元胞）在各个方向上重复叠合而成。但在模拟中只需保留基本单 元，所有其它单元与基本单元由平移对称性关联。
经典的分子动力学方法是 Alder 等于 1957 年提出并首先在“硬球”液体模型下 应用，发现了由 Kirkwood 在 1939 年根据统计力学预言的“刚性球组成的集合系统 会发生有液相到结晶相的转变”。后来人们称这种相变为 Alder 相变。Rahman 于 1963 年采用连续势模型研究了液体的分子动力学模拟。1972 年 Less 等发展了该方法并 扩展了存在速度梯度的非平衡系统。1980 年 Andersen 等创造了恒压分子动力学方 法。1983 年 Gillan 等将该方法推广到具有温度梯度的非平衡系统，从而形成了非 平衡系统分子动力学方法体系。1984 年 Nose 等完成了恒温分子动力学方法的创建。 1985 年针对势函数模型化比较困难的半导体和金属等，Car 等提出了将电子论与分 子动力学方法有机统一起来的第一性原理分子动力学方法。1991 年 Cagin 等进一步 提出了应用于处理吸附问题的巨正则系综分子动力学方法。20 世纪 80 年代后期， 计算机技术飞速发展，加上多体势函数的提出与发展，使分子动力学模拟技术有了 进一步的发展。 1、分子动力学的运动方程：
分子动力学模拟就是用计算机方法来表示统计力学，作为实验的一个辅助手段。 分子模拟就是对于原子核和电子所构成的多体系统，求解运动方程(如牛顿方程、哈
密顿方程或拉格朗日方程)，其中每一个原子核被视为在全部其它原子核和电子作用 下运动，通过分析系统中各粒子的受力情况，用经典或量子的方法求解系统中各粒 子在某时刻的位置和速度，以确定粒子的运动状态，进而计算系统的结构和性质。 该模拟技术主要涉及粒子运动的动力学问题,与蒙特卡罗模拟方法(简称 MC)相比, 分子动力学是一种“确定性方法”,它所计算的是时间平均,而 MC 进行的是系综平均。 然而按照统计力学各态历经假设,时间平均等价于系综平均。因此,两种方法严格的 比较计算能给出几乎相同的结果。
分子动力学模拟方法的基本原理与应用
摘要: 介绍了分子动力学模拟的基本原理及常用的原子间相互作用势, 如 Lennard-Jones 势; 论述 了几种常用的有限差分算法, 如 Verlet 算法; 说明了分子动力学模拟的几种系综及感兴趣的宏观 统计量的提取。 关键词: 分子动力学模拟; 原子间相互作用势; 有限差分算法;
2、运动方程的数值积分： 计算机模拟方法的基点是利用现代计算机高速和精确的优点 , 对几百个以至
上千个分子的运动方程进行数值积分 有许多不同的积分方法 , 它们的效率和方便 程度各异 问题基本上就是用有限差分法来对二阶常微分方程进行积分常用的有以 下几种方法： 、Verlet 算法［1］：
Verlet 算法是在 60 年代后期出现的，对扩散分子的质心运动的积分是最稳定 的也是最常用的数值方法。它运用 t 时刻的位置和加速度以及 t 时刻的位置来预测 t＋δt 位置，其积分方案，以三阶 Taylor 展开为基础，由以下方程给出：
r(t＋δt)＝2r(t)-r(t-δt)+δt2α(t) 这里，为简单计，省略了i。速度可按微分的基本法则得出： V(t)=[r(t+δt)-r(t-δt)]/2δt。这种算法的优点是占有计算机的内存小，并且 很容易编程。但它的缺点是位置r(t+δt)要通过小项(δt2)与非常大的两项2r(t)和 r(t-δt)的差相加得到，这容易造成精度损失。并且从式中可以看出，这种算法不 是一个自启动算法，新位置必须由t和t-δt时刻的位置得到。
多体势是在 20 世纪 80 年代初期开始出现的，1984 年 Daw 和 Baskes 首次提出 了原子嵌入（EAM）势[4]。此势的基本思想是：把晶体的总势能分成位于晶格点阵上 的原子核之间的相互作用对势和原子核镶嵌在电子云背景中的嵌入能（多体相互作 用）两部分，其中，对势和多体作用势的函数形式往往根据经验选取。基于 EAM 势 的势函数还有很多种，这些多体势大多用于金属的微观模拟，此外，还有许多形式 的多体势函数形式，1987 年 Jacobsen 等人在等效介质原理的基础上提出的另一种 多体势函数形式，由于其简单、有效，也得到了广泛地应用。 5、分子动力学模拟的系综：
、对势：
在分子动力学模拟的初期，经常采用的就是对势。对势可以分为间断对势[1]和 连续对势，而连续对势主要有以下几种：Lennard-Jones（L-J）势、Born-lande（B-L） 势、Morse 势和 Johnson 势等，其中，L-J 势是为描述惰性气体分子之间的相互作用 而建立的，因此它表达的力比较弱，描述材料的行为也比较柔韧，也可以用来描述 过渡金属原子之间的相互作用。B-L 势是用来描述离子晶体离子之间的相互作用的。 Morse 势和 Johnson 势多用于描述金属原子之间的相互作用。对势虽然简单，得到 的结果往往也符合某些宏观物理规律，但其缺点是必然导致 Cauchy 关系，所以，对 势实际上不能准确地描述晶体的弹性性质。 、多体势：
分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子 模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的 不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果 为基础进一步计算体系的热力学量和其他宏观性质。
从 统 计 物 理 学 中 衍 生 出 来 的 分 子 动 力 学 模 拟 方 法 (Molecular Dynamics Simulation, MDS) , 实践证明是一种描述纳米科技研究对象的有效方法, 得到越来 越广泛的重视。所谓分子动力学模拟, 是指对于原子核和电子所构成的多体系统, 用计算机模拟原子核的运动过程, 从而计算系统的结构和性质, 其中每一个原子核 被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动。它被认 为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段, 称之为“计算 机实验”手段, 在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。
ai=-(dV/dr+mivi )/mi，d
。式中，g 为体系的自由度，Q
为一个可调参量，表征着假想项的质量，T 为温度。 在这三种控温方法中，速度标定是一种非真实的物理效应。但这种方法可以使
系统很快达到平衡，在经典分子动力学方法中，这是一种比较常用的控温方法。 Berendsen 控温方法是通过系统和恒温槽进行热交换来控制温度，此方法的优点是 非常简单，应用起来非常方便。Nose-Hoover 控温方法是基于统计力学而提出来的 一种控温方法，当系统和恒温槽进行热交换时，在系统中粒子出现的几率遵从统计 力学规律，所以这是一种真实的物理效应。 、控压方法:
（3）Nose-Hoover 方法：这种方法是通过改变模拟体系的 Hamiltonian 来实现
控温的，因而有更强的物理意义，其基本思想就是在 Hamiltonian 加入一个假想的
项来代表一个恒温源，具体做法如下：
，
式中，S 和 分别是假想项的坐标和动量。这样体系的微分方程就变为：vi=dri/dr，
、Gear的预侧-校正算法：
这种算法分为三步来完成：首先，根据 Taylor 展开，预测新的位置、速度和加 速度。然后，根据新的计算的力计算加速度。这个加速度再由与 Taylor 级数展开式 中的加速度进行比较，两者之差在校正步里用来校正位置和速度项。这种方法的缺 点就是占有计算机的内存大。 、“蛙跳”（Leap-frog）算法[2]：
温度的变化为ΔT=(λ2-1)T(t)，其中，
，T 为所控制体系的温度。
（2）Berendsen 恒温槽方法：这种控温方法假设所模拟的体系与一个恒温槽连
载一起，则两者之间就可以通过热交换而使模拟体系达到恒温的目的，方法如下，
定义一个参数 λ：
，式中τ表征系Leabharlann 与恒温槽之间的热交换速率，Δt 为 MD 的时间步长，那么通过 vnew=λvold 校正即可保持体系的温度在 T0 附近震动，而参数 （通常取为~）则可用于控制这个震动幅度。
平衡态分子动力学模拟是在一定的系综瞎进行的，经常用的平衡系综是 NVT 或 NPT 系综。在这两种系综中，牵涉到控制温度和压力的几种技术，分别如下： 、控温方法：
（1）速度标定：系统的温度与动能存在如下关系：
，
式中，N 是原子数，Nc 是约束数，KB 是 Boltzmann 常数，vi 是原子 i 的速度。由于系 统的温度和动能存在这样的关系，所以一种最简单和最直观的方法是直接对速度进 行标定。这种方法的基本思想是：如果 t 时刻的温度是 T(t)，速度乘以因子λ后，
温方法类似，如下：
，式中 是一个可调参数，
在处理粒子之间的相互作用时，通常采用“最小影像”约定。这个约定是在由 无穷重复的 MD 基本模拟计算元胞中，一个粒子只与它所在的基本元胞内的另外 N-1 个（设在此元胞内有 N 个粒子） 中的每个粒子或其最邻近影像粒子发生相互作用。 实际上，这个约定就是通过满足不等式条件 rc<L/2 来截断位势（rc 为截断半径，L 是元胞的边长）通常 L 的数值应当选得很大，以避免有限尺寸效应，但这样会增大 计算量，同城采用对相互作用势的修正来近似处理。 4、势函数：
科学工作者在长期的科学研究实践中发现，当实验研究方法不能满足研究工作 的需求时，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；尽 管计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从 而促进理论和实验的发展。特别是在材料形成过程中许多与原子有关的微观细节， 在实验中基本上是无法获得的，而在计算机模拟中即可以方便地得到。这种优点使 分子动力学模拟在材料研究中显得非常有吸引力。
用 MD 方法研究压力的诱导相变和结构重构，在等压下模拟比在等体积下更容易 实现。在等压模拟下，可以通过改变模拟元胞的三个方向或一个方向的尺寸来实现 体积的变化，类似于温度控制的方法，也有许多方法用于压力控制，在实际应用中， 有以下两种较为常用的压力控制方法。
（1） Berendsen 方法：这种控压方法的基本思路与基本控制方程均与他的控
分子动力学模拟的出发点是假定粒子的运动可以用经典动力学来处理，对一个 由Ｎ个粒子构成的孤立体系，粒子的运动由牛顿运动方程决定，也就是： mid2ri/dt2=-▽iV(r1,r2,rN),式中，mi，ri 分别为第 i 个原子的质量和位置。▽ i=-?/?ri，V(r1,r2,rN)为体系所处的势。
Hockey 提出的 Leap-frog 算法是 Verlet 算法的变化，这种方法设计半时间间 隔的速度，即：r(t+δt)=r(t)+ δtv(t+δt/2)，v(t+δt/2)=v(t-δt/2)+ δtα(t)。
t 时刻的速度由下式给出：v(t)=[v(t+δt/2)+v(t-δt/2)]/2。这种算法与 Verlet 算法相比有两个优点：（1）包括显速度项；（2）收敛速度快，计算量小。这种算法 明显的缺陷是位置和速度不同步。
除了上述提及的几种方法外, 还有 Beeman 算法、Rahman 等。 3、周期性边界条件和长程力：
即使是使用现代的巨型计算机，MD 方法还是只能用于粒子数大约是几百到几千 的系统。这就引起一个问题：用这样少量的粒子，如何来模拟宏观体系为了解决这 个问题，引入了周期性边界条件[3]。采用这种方法，模拟体系实际上是由基本单元 （也称为模拟计算元胞）在各个方向上重复叠合而成。但在模拟中只需保留基本单 元，所有其它单元与基本单元由平移对称性关联。
经典的分子动力学方法是 Alder 等于 1957 年提出并首先在“硬球”液体模型下 应用，发现了由 Kirkwood 在 1939 年根据统计力学预言的“刚性球组成的集合系统 会发生有液相到结晶相的转变”。后来人们称这种相变为 Alder 相变。Rahman 于 1963 年采用连续势模型研究了液体的分子动力学模拟。1972 年 Less 等发展了该方法并 扩展了存在速度梯度的非平衡系统。1980 年 Andersen 等创造了恒压分子动力学方 法。1983 年 Gillan 等将该方法推广到具有温度梯度的非平衡系统，从而形成了非 平衡系统分子动力学方法体系。1984 年 Nose 等完成了恒温分子动力学方法的创建。 1985 年针对势函数模型化比较困难的半导体和金属等，Car 等提出了将电子论与分 子动力学方法有机统一起来的第一性原理分子动力学方法。1991 年 Cagin 等进一步 提出了应用于处理吸附问题的巨正则系综分子动力学方法。20 世纪 80 年代后期， 计算机技术飞速发展，加上多体势函数的提出与发展，使分子动力学模拟技术有了 进一步的发展。 1、分子动力学的运动方程：
分子动力学模拟就是用计算机方法来表示统计力学，作为实验的一个辅助手段。 分子模拟就是对于原子核和电子所构成的多体系统，求解运动方程(如牛顿方程、哈
密顿方程或拉格朗日方程)，其中每一个原子核被视为在全部其它原子核和电子作用 下运动，通过分析系统中各粒子的受力情况，用经典或量子的方法求解系统中各粒 子在某时刻的位置和速度，以确定粒子的运动状态，进而计算系统的结构和性质。 该模拟技术主要涉及粒子运动的动力学问题,与蒙特卡罗模拟方法(简称 MC)相比, 分子动力学是一种“确定性方法”,它所计算的是时间平均,而 MC 进行的是系综平均。 然而按照统计力学各态历经假设,时间平均等价于系综平均。因此,两种方法严格的 比较计算能给出几乎相同的结果。
分子动力学模拟方法的基本原理与应用
摘要: 介绍了分子动力学模拟的基本原理及常用的原子间相互作用势, 如 Lennard-Jones 势; 论述 了几种常用的有限差分算法, 如 Verlet 算法; 说明了分子动力学模拟的几种系综及感兴趣的宏观 统计量的提取。 关键词: 分子动力学模拟; 原子间相互作用势; 有限差分算法;
2、运动方程的数值积分： 计算机模拟方法的基点是利用现代计算机高速和精确的优点 , 对几百个以至
上千个分子的运动方程进行数值积分 有许多不同的积分方法 , 它们的效率和方便 程度各异 问题基本上就是用有限差分法来对二阶常微分方程进行积分常用的有以 下几种方法： 、Verlet 算法［1］：
Verlet 算法是在 60 年代后期出现的，对扩散分子的质心运动的积分是最稳定 的也是最常用的数值方法。它运用 t 时刻的位置和加速度以及 t 时刻的位置来预测 t＋δt 位置，其积分方案，以三阶 Taylor 展开为基础，由以下方程给出：
r(t＋δt)＝2r(t)-r(t-δt)+δt2α(t) 这里，为简单计，省略了i。速度可按微分的基本法则得出： V(t)=[r(t+δt)-r(t-δt)]/2δt。这种算法的优点是占有计算机的内存小，并且 很容易编程。但它的缺点是位置r(t+δt)要通过小项(δt2)与非常大的两项2r(t)和 r(t-δt)的差相加得到，这容易造成精度损失。并且从式中可以看出，这种算法不 是一个自启动算法，新位置必须由t和t-δt时刻的位置得到。
多体势是在 20 世纪 80 年代初期开始出现的，1984 年 Daw 和 Baskes 首次提出 了原子嵌入（EAM）势[4]。此势的基本思想是：把晶体的总势能分成位于晶格点阵上 的原子核之间的相互作用对势和原子核镶嵌在电子云背景中的嵌入能（多体相互作 用）两部分，其中，对势和多体作用势的函数形式往往根据经验选取。基于 EAM 势 的势函数还有很多种，这些多体势大多用于金属的微观模拟，此外，还有许多形式 的多体势函数形式，1987 年 Jacobsen 等人在等效介质原理的基础上提出的另一种 多体势函数形式，由于其简单、有效，也得到了广泛地应用。 5、分子动力学模拟的系综：
、对势：
在分子动力学模拟的初期，经常采用的就是对势。对势可以分为间断对势[1]和 连续对势，而连续对势主要有以下几种：Lennard-Jones（L-J）势、Born-lande（B-L） 势、Morse 势和 Johnson 势等，其中，L-J 势是为描述惰性气体分子之间的相互作用 而建立的，因此它表达的力比较弱，描述材料的行为也比较柔韧，也可以用来描述 过渡金属原子之间的相互作用。B-L 势是用来描述离子晶体离子之间的相互作用的。 Morse 势和 Johnson 势多用于描述金属原子之间的相互作用。对势虽然简单，得到 的结果往往也符合某些宏观物理规律，但其缺点是必然导致 Cauchy 关系，所以，对 势实际上不能准确地描述晶体的弹性性质。 、多体势：
分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子 模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的 不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果 为基础进一步计算体系的热力学量和其他宏观性质。
从 统 计 物 理 学 中 衍 生 出 来 的 分 子 动 力 学 模 拟 方 法 (Molecular Dynamics Simulation, MDS) , 实践证明是一种描述纳米科技研究对象的有效方法, 得到越来 越广泛的重视。所谓分子动力学模拟, 是指对于原子核和电子所构成的多体系统, 用计算机模拟原子核的运动过程, 从而计算系统的结构和性质, 其中每一个原子核 被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动。它被认 为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段, 称之为“计算 机实验”手段, 在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。
ai=-(dV/dr+mivi )/mi，d
。式中，g 为体系的自由度，Q
为一个可调参量，表征着假想项的质量，T 为温度。 在这三种控温方法中，速度标定是一种非真实的物理效应。但这种方法可以使
系统很快达到平衡，在经典分子动力学方法中，这是一种比较常用的控温方法。 Berendsen 控温方法是通过系统和恒温槽进行热交换来控制温度，此方法的优点是 非常简单，应用起来非常方便。Nose-Hoover 控温方法是基于统计力学而提出来的 一种控温方法，当系统和恒温槽进行热交换时，在系统中粒子出现的几率遵从统计 力学规律，所以这是一种真实的物理效应。 、控压方法:
（3）Nose-Hoover 方法：这种方法是通过改变模拟体系的 Hamiltonian 来实现
控温的，因而有更强的物理意义，其基本思想就是在 Hamiltonian 加入一个假想的
项来代表一个恒温源，具体做法如下：
，
式中，S 和 分别是假想项的坐标和动量。这样体系的微分方程就变为：vi=dri/dr，
、Gear的预侧-校正算法：
这种算法分为三步来完成：首先，根据 Taylor 展开，预测新的位置、速度和加 速度。然后，根据新的计算的力计算加速度。这个加速度再由与 Taylor 级数展开式 中的加速度进行比较，两者之差在校正步里用来校正位置和速度项。这种方法的缺 点就是占有计算机的内存大。 、“蛙跳”（Leap-frog）算法[2]：