硅谷网：探秘分子动力学模拟方法的优势和发展方向 |
| 2013-05-20 18:06 作者：周晓平 田壮壮 忽晓 来源：《硅谷》杂志 HV： 编辑： 【搜索试试】
|
|
硅谷网文 据《硅谷》杂志2012年第23期刊文称，主要介绍了分子动力学模拟的基本原理，阐述分子动力学方法的运动方程、数值解法、势函数、边界条件、适用系综以及体系相关性质的计算。最后指出分子动力学模拟方法的优势和发展方向。
1分子动力模拟计算的基本原理
分子动力计算的基本原理，即为利用牛顿运动定律。在分子动力模拟中，体系原子的一系列位移是通过对牛顿运动方程的积分得到的，结果是一条运动轨迹，它表明了系统内原子的位置与速度是如何随时间而发生变化。
先由系统中各分子位置计算系统的势能，按照经典力学，系统中任一原子i所受的力为势能的梯度：
（1）
由牛顿第二定律可得i原子的加速度为：
（2）
将牛顿运动定律方程式对时间积分，可预测i原子经过时间t后的速度与位置：
，
，
（3）
式中，及分别是粒子i的位置与速度，上标“0”为各物理量的初始值[]。
2牛顿运动方程的数值解法
为了得到原子的运动轨迹，必须解式（3）的牛顿运动方程，可采用有限差分法。有限差分法的基本思想就是将积分分成很多小步，每一小步的时间固定为。常用的有以下几种算法：①Verlet算法；②Velocity-Verlet算法；③leap-frog算法（蛙跳算法）；④Beeman算法；⑤Gear算法。leap-frog算法和Gear算法由于使用简便，准确性及稳定性高，节省储存空间等作者：photon优点，已被广泛采用。
2.1leap-frog算法
leap-frog算法速度与位置的数学式为：
（4）
为了执行leap-frog算法，必须首先由t-0.5时刻的速度与t时刻的加速度计算出，然后由方程
（5）
计算出位置。时间为t时的速度可由式（6）算出，即：
（6）
利用蛙跳法计算仅需储存两种资料，可节省储存空间。其缺点是位置与速度不同步。这意味着在位置一定时，不可能同时计算动能对总动能的贡献。
2.2Gear算法[1]
Gear所提出的一种利用数值解的方法，称为校正预测法（predictor-correctormethod）。时间t+时的位置、速度等可由时间t的泰勒展开式预测得到：
（7）
式中的的1次、2次、3次微分。式（7）所产生的速度、加速度不是由牛顿运动方程解得的，所以并非完全正确。可由所预测的位置计算所受的力及正确的加速度。设正确的加速度与预测的加速度之间的误差为：
（8）
可得各量的校正式为：
（9）
式中，均为常数。以上为Gear的一次预测校正法，也可将此计算推展至更高次的校正。
3势函数
势函数表明了原子间的相互作用。针对不同的计算物质，不同的模拟目的，势函数有不同的形式。计算结果的可靠性与势函数密切相关。在分子动力学发展初期，主要采用对势。随着模拟体系的复杂性，逐渐出现了多体势，以弥补对势的不足。
3.1对势
主要是Lennard-Jones势（L-J势），又叫12－6势能，它的数学表达式是：
（10）
式中，r为原子对间的距离，、是势能参数。在L-J势能中，项是排斥项，项是吸引项。当r很大时，L-J势能趋近于零，表示当原子对相距很远时，彼此之间已经没有非键结作用了。这个形式的势函数表达的作用力较弱，适合描述惰性气体的固体和液体[]。
3.2多体势
Daw和Baskes首次提出了嵌入原子势（EAM），其体系的总势能可以表示为：
（11）
式中，第一项是嵌入能，表示原子核镶嵌在电子云背景中的嵌入能；第二项是对势能，表示位于晶格点阵上的原子核之间的相互作用。多体势大都用来模拟金属的微观性质[2]。后来Baskes等人修正了嵌入原子法，提出MEAM势，可以描述共价键材料。
4周期性边界条件与最近镜像
执行分子动力计算时，通常选取一定数目N的分子或原子，将其放在一个立方体的盒子中，该盒子即为模拟系统。而计算机最多只能模拟几百到几千个粒子的系统，为了能够用较少的粒子数目来模拟真实的宏观体系，引入三维周期性边界条件。在三维体系中，每一个单胞会被其它26个单胞所包围。当模拟的单胞中一个粒子由于力的作用离开这个单胞的时候，就会有另一个和它对应的粒子运动到这个单胞中来，这样，模拟的整个体系的粒子数就保持不变，密度也不变，符合实际的要求[]。
计算系统中分子间的作用力时，采取最近镜像的方法。计算分子A与B的作用力，是取与分子A和其最近的距离镜像分子B，而非计算系统中的分子A与分子B。因为在计算中利用了最近镜像的概念，因此必须采用截断半径的方法计算非键结的远程作用力。定义为截断半径，则当时，势能值趋近于零，表示分子间的范德瓦耳斯作用力可忽略不计。截断半径最大不能超过盒子边长的一半，即。我们在进行模拟计算时，一般的原子所选取的半径为10～15Å[1]。
5各种系综的分子动力学模拟
采用分子动力学方法模拟，必须在一定的系综下进行。经常用到的系综有正则系综（NVT）、微正则系综（NVE）、等温等压系综（NPT）、等压等焓系综（NPH）等。
5.1正则系综
正则系综是针对一个粒子数N、体积V、温度T保持不变，同时总动量P为零的系综。为了保持系统温度不变，通常运用的方法是让系统与外界的热浴处于热平衡状态。由于温度与系统的动能有直接的关系，正则系综的分子动力学模拟就是在运动方程上加上动能恒定的约束。
5.2微正则系综
微正则系综中，粒子数N、体积V、能量E都保持不变，它是一个孤立的、保守的系统的系综，在这种系综中，系统沿着相空间中的恒定能量轨道演化。
5.3等温等压系综
等温等压系综是常见的系综，即系统的粒子数N、温度T、压强P保持不变。由于系统的压力P与体积V是共轭量，对压力的调节可以通过改变系统的体积来实现。
5.4等压等焓系综
等压等焓系综中，系统粒子数N、压强P、焓值H保持不变，焓值通过H=E+P得到。在此系综下进行模拟时，一定要保持压力和焓值为固定值。
6物理量的提取
6.1径向分布函数
径向分布函数（radialdistributionfunction，RDF）是表征系统结构的特征量，尤其是对液体。径向分布函数g（r）定义为：
（12）
式中，为系统密度，dN为距离参考分子中心由球壳内的平均粒子数。
径向分布函数可以理解为所模拟体系的区域密度与平均密度的比[1]。当r值小时，距离参考分子很近，区域密度不同于平均密度；当r值大时，距离参考分子很远，区域密度与系统的平均密度相同，径向分布函数接近于1。
6.2配位数
配位数（coordinationnumber）是指某个原子的最近邻原子个数，它显示了该原子周围的原子分布的密度大小，是微观结构分析的辅助物理量[]。其定义为：
（13）
式中，是对应点的数密度，积分的上下限由径向分布函数g（r）得到，当时，得到第一配位圈的配位数。
6.3扩散系数
分子动力学模拟中，采用两种方法来计算体系分子的扩散系数D，分别是Einstein法和Green-Kubo法。Eintein法是利用对均方位移求斜率来求扩散系数，计算公式为：
（14）
Green-Kubo法是通过对速度自相关函数（）的积分获得扩散系数，计算公式为：
（15）
其中，表示粒子i的扩散系数，、分别表示粒子在时刻和0时刻的位置矢量。
7结束语
分子动力学模拟方法不仅可以提供物质体系的宏观和微观性质，还可以模拟极端条件（如高温、高压[]）下物质的各种性质，完善实验数据，为广大研究者提供重要参考。随着计算机技术的长足发展，利用计算机对物质体系进行模拟得到了广泛应用。因为势函数对模拟结果的精度影响较大，提高势函数的精确性会是分子动力学模拟未来的发展重点。
基金项目：河南省教育厅科学技术研究重点项目（12B510029）
作者简介：
周晓平（1984-），汉族，女，硕士，主要从事原子与分子团簇方面的研究。 |
|
|
|
【对“硅谷网：探秘分子动力学模拟方法的优势和发展方向”发布评论】 |
版权及免责声明：
① 本网站部分投稿来源于“网友”，涉及投资、理财、消费等内容，请亲们反复甄别，切勿轻信。本网站部分由赞助商提供的内容属于【广告】性质，仅供阅读，不构成具体实施建议，请谨慎对待。据此操作，风险自担。
② 内容来源注明“硅谷网”及其相关称谓的文字、图片和音视频，版权均属本网站所有，任何媒体、网站或个人需经本网站许可方可复制或转载，并在使用时必须注明来源【硅谷网】或对应来源，违者本网站将依法追究责任。
③ 注明来源为各大报纸、杂志、网站及其他媒体的文章，文章原作者享有著作权，本网站转载其他媒体稿件是为传播更多的信息，并不代表赞同其观点和对其真实性负责，本网站不承担此类稿件侵权行为的连带责任。
④ 本网站不对非自身发布内容的真实性、合法性、准确性作担保。若硅谷网因为自身和转载内容，涉及到侵权、违法等问题，请有关单位或个人速与本网站取得联系（联系电话：01057255600），我们将第一时间核实处理。
|
|
|
|
