分子动力学模拟的计算及应用

摘 要： 主要介绍分子动力学模拟的来源，从分子动力学的运动方程出发，详细介绍分子动力学模拟过程中的数据分析方法：有限差分法，分子动力学模拟所适用的系综，以及最后对于模拟数据的热力学性质的计算。

关键词： 分子动力学模拟；有限差分法；热力学性质

中图分类号：O561 文献标识码：A 文章编号：1671－7597（2012）0220149－01

所谓分子动力学模拟是利用计算机技术来模拟或仿真已有的理论知识，进一步来分析微观方面分子或原子的运动。随着人们对牛顿力学的不断研究，进一步从原子物理和统计物理的角度来从微观方面分析宏观的问题。目前所具有的配置低的计算机硬件无法满足后来学者的需求，此时人们就进一步的应用国外先进的配置高的计算机来详细及的对凝聚态物质分子方面微观结构的模拟计算。与此同时还研究出了两种用数值方法运算统计求和的方法：分子动力学（MD）方法和蒙特卡罗（MC）方法。MD和MC都是从热力学的角度来研究凝聚态物质的分子或原子的运动行为。但MD还不同于MC，MD是以组成系统的大量的原子或分子来作为研究对象研究的。1957年，Alder和Wainwright两人首先在研究硬球模型下的液体模型时首次提出的。他们发表的第一篇有关分子动力学的文章，里面主要研究了硬球体系的相位图。然后通过模拟得到的数据来从多个方面研究这一系统中微观物质的宏观性质。这些的到的模拟数据还可以通过作图来准确的描述整个模拟过程中的状态以及系统中分子的性质。结果与理论结果相近似或相吻合，有利于对新的概念和理论的提取。目前，计算机所做的分子动力学模拟已成为物理学、化学、材料科学及制药方面研究必不可少的工具，并且分子动力学模拟（MD模拟）以应用于模拟原子的扩散、相变、薄膜生长和表面缺陷等过程。

1 分子动力学的运动方程

分子动力学研究的对象主要是由一定量的分子组成的系综，假设现在研究的分子数为N，则牛顿运动方程来描述分子的运动，也就是：

式中mi与ri分别为第i个粒子的质量和位置。体系的所处的势为

2 有限差分算法

2.1 维里算法Verlet[7]

位置：r（t+Δt）=2r（t）-r（t-Δt）+（Δt）2a（t）

用一个简单的方法将t+Δt和t-Δt时刻粒子的位置坐标分别做泰勒展

其中公式中m表示院子的质量，V（t）表示院子的速度，F（t）为原子所受到的力，而a（t）表示原子的加速度。

2.2 Velocity-Verlet算法[1]

Velocity-Verlet算法可以同时给出位置、速度与加速度，并且精度非常高，而且计算量不是太大，对时间是可逆的。表示为：

2.3 Gear算法：（校正预测算法）

MD模拟中，是在微正则系综中，想要能量保持守恒就尽可能用较大得时间不长来模拟，然后再应用矫正预测法，实验结果更准确。这种算法从Taylor展开开始，用{rin}得到{rin+i}来预测新的位置和速度。进一步得到加速度a（t+δt），然后再让次加速度与Taylor级数展开式中的加速度ac（t+δt）进行比较。这样得到的结果更精确、更可靠。

2.4 Leap-frog（蛙跳算法）

维里算法速度项中含有1/Δt项，由于实际计算中一般选取的Δt都很

。T时刻的速度由下式给出：

。此算法与与Verlet算法相比有两个优点：1）包括显速度项；2）收敛速度快，计算量小。

3 分子动力学模拟的不同系综

3.1 微正则系综（NVE）

微正则系综是孤立的、保守的系统，整个系综与外界既无粒子交换又无能量交换，即所模拟系统中的粒子数N、体积V、能量E都保持不变。而且由于整个系统没有发生运动，故整个系统的总动量为零。

3.2 正则系综（NVT）[8]

在正则系综中，系统本身的粒子数N、体积V和温度T都始终保持不变。在恒温的情况下，系统的总能量并不是个守恒量，系统要靠与外界发生能量交换来保持系统整体的温度不变，一般采用的方法是让整个系统与外界的恒温大热源相互接触，进行传热。

3.3 等温等压系综（NPT）

分子动力学模拟中的等温等压系综，它的粒子数N、温度T和压力P都保持恒定不变。温度控制和上面一样，通过与恒温大热源接触来保持恒温。

3.4 等焓等压系综（NHP）

所谓等焓等压系综就是系统的粒子数N、压力P和焓值H都保持不变。模拟时要保持压力与焓值为固定值。这种系综在实际的分子动力学模拟中很少见。

4 热力学性质的计算

1）径向分布函数（RDF）：径向分布函数表示粒子的聚集状况，在研究液态和非晶态的结构中常用。表达式为：

几率，p是系统的平均数密度，R是原子位置，δ是Dirac符号，N为原子数。

2）静态结构因子（SSF）：是判断分子无序化程度的物理量，其表达

的位置矢量。理想晶体的静态结构因子为1，理想流体怎为0。

中：ri（O）为原子i在零时刻的位移；ri（t）为原子i在t时刻的位移。由爱因斯坦扩散定律，均方位移随时间的变化表征了液态金属原子的扩散行为。它与扩散系数（D）存在如下关系： 式中C为常数。

3）局部晶序分析（LCO）：又称对分析技术法，用来计算不同温度下原子间键合类型及指数。根据这种方法，用4位数ijkl描述原子所属的状态：i代表两个原子的成键关系，i-1为成键，i-2为未成键；j代表成键两原子的共用最近邻原子数；k代表共用最近邻原子之间的成键数；对前三

数相同而结构不同的用不同的l值来区分。

4）配位数（CN）：配位数作为一种分析结构的辅助手段，是指某原子的第一最近邻原子的个数。进一步来说明某一原子周围的原子分布的密度。

参考文献：

[1]杨萍、孙益民，分子动力学模拟方法及其应用，安徽师范大学学报（自然科学版），2009，1，32（1）.

[2]崔守鑫、胡海泉、肖效光，分子动力学模拟基本原理和主要技术，聊城大学学报（自然科学版），2005，3，15（1）.

[3]文玉华、朱如曾、周富信等，分子动力学模拟的主要技术，力学进展，2003，33（1）.

[4]樊康旗、贾建授，经典分子动力学模拟的主要技术，器件与技术.

[5]申海兰、赵靖松，分子动力学模拟方法概述，装备制造技术，2007（10）.

[6]赵素、李金富等，分子动力学模拟及其在材料科学中的应用，材料导报，2007，4，21（4）.