分子动力学(Molecular Dynamics，简称MD)模拟就是在给定的原子间相互作用势的情况下，结合牛顿力学的基本原理，通过对原子核和电子所构成的多体系统求解牛顿运动方程得到体系中原子核的运动过程，进而计算出体系的结构和性质，其中每一个原子核可视为在全部其他原子核和电子的作用下运动. 1957年Alder和Wainwright^[1] 通过计算机模拟的方法，研究从32个到500个刚性小球的分子系统的运动，开创利用MD模拟方法研究物质宏观性质的先例. Lees等^[2]在1972年和Gillan^[3]在1983年先后将该方法扩展到存在速度和温度梯度的非平衡体系；Andersen等^[4]在1980年和Nosé 等^[5]在1984年又分别完成恒压分子动力学和恒温分子动力学方法的创建；1985年Hoover等^[6]对Nosé的方法进行改进使得该方法更加稳定可靠；同年，Car等^[7]针对势函数模型化比较困难的半导体和金属等提出将电子论和分子动力学方法结合的第一性原理分子动力学方法；1991年，Cagin等^[8]进一步提出应用于处理吸附问题的巨正则系统(体系的体积、温度、化学势能为定值)的MD方法. 20世纪80年代后期，计算机的飞速发展和多体势函数的提出与应用，使得MD模拟方法更加成熟，由于很多实验的特殊环境需求(超临界、深过冷、高温等)，常规实验难以满足，越来越多的MD模拟软件(Moldy、LAMMPS、DL_POLY、Materials Studio)被开发出来，而且通过这些软件的模拟，问题得到很好地解决.

1 分子动力学模拟的原理

分子动力学模拟是在原子及分子水平上对多体问题进行求解的计算机模拟方法.在模拟计算中根据2个基本假设^[9]：①体系中所有粒子的运动都遵循牛顿运动定律；②粒子之间的相互作用满足叠加定理.这2条基本假设忽略量子效应和多体作用，与真实的物理系统存在一定的差别，故分子动力学模拟属于近似计算.MD模拟的出发点是假设N为模拟系统的原子数，m_i为第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₁，r₂，k_N)=$\sum\limits_{i=1}^{N}{{}}\sum\limits_{j=1}^{N}{{}}$Φ(r_ij)，所给定物理量都是随时间而变化的，即A=A(t)，控制方程如下：

${{m}_{i}}{{r}_{i}}={{F}_{i}}=\sum\limits_{j\ne 1}^{{}}{{}}{{f}_{ij}}$

(1)

${{F}_{i}}=-{{\Delta }_{{{r}_{i}}}}U\left( {{r}_{1}},{{r}_{2}},{{k}_{{{r}_{N}}}} \right)$

(2)

由式(1)和式(2)共同建立一个封闭的线性微分方程组，并且通过给定的初始位置和速度，可以求得系统中所有原子任意时刻的位置r_i(t)和速度V_i(t).

2 运动方程的数值积分

MD模拟计算中的大多数算法都是有限差分法对牛顿运动方程进行积分的.常用的有以下几种.

2.1 Verlet算法

Verlet算法^[10]是在20世纪60年代后期提出的，由于它的计算过程比较简单，对扩散分子的质心运动的积分也最稳定，所以在MD模拟中运用最为广泛.它运用原子在t时刻的位置r(t)和加速度a(t)以及t-δ(t)时刻的位置r(t-δt),计算出原子在t+δ(t)时刻的位置r(t＋δt),推导方程如下：

$r\left( t+\delta t \right)=r\left( t \right)+\delta tv\left( t \right)+\frac{1}{2}\delta {{t}^{2}}a\left( t \right)+\ldots $

(3)

$r\left( t-\delta t \right)=r\left( t \right)-\delta tv\left( t \right)+\frac{1}{2}\delta {{t}^{2}}a\left( t \right)-\ldots $

(4)

式(3)和式(4)相加得：

$r\left( t+\delta t \right)=2r\left( t \right)-r\left( t-\delta t \right)+\delta {{t}^{2}}a\left( t \right)+O\left( \delta {{t}^{4}} \right)$

式(3)和式(4)相减得：

$v\left( t \right)=\frac{r\left( t+\delta t \right)-r\left( t-\delta t \right)}{2\delta t}+O\left( \delta {{t}^{2}} \right)$

这种算法的优点是占用计算机的内存小，包括两组位置坐标r(t),r(t-δt)和一组加速a(t).但它的缺点是位置r(t＋δt)需要通过(δt²)与2r(t)和r(t-δt)两项的差相加得到，降低计算精度；其速度计算为二级近视，坐标却达四级近视，相差两个近视等级；另外一个是其坐标和速度的计算数据没有关联，每一时间的速度计算都滞后于坐标计算，且精度比较高.

2.2 速度Verlet算法（velocity Verlet algorithm)

Swope提出的速度Verlet算法^[11]可以同时给出位置、速度与加速度，且不消耗精度，是Verlet算法的改进形式.其推导过程如下： 已知

$v\left( t+\delta t \right)=v\left( t \right)+\delta ta\left( t \right)+\frac{1}{2}\delta {{t}^{2}}b\left( t \right)+O\left( \delta {{t}^{3}} \right)$

(5)

又

$a\left( t+\delta t \right)=a\left( t \right)+\delta tb\left( t \right)+\ldots ,=>\delta tb\left( t \right)=a\left( t+\delta t \right)-a\left( t \right)$

(6)

将式(6)代入式(5)可得：

$v\left( t+\delta t \right)=v\left( t \right)+\delta ta\left( t \right)+\frac{\delta t}{2}\left[ a\left( t+\delta t \right)-a\left( t \right) \right]+O\left( \delta {{t}^{3}} \right)$

即

$v\left( t+\delta t \right)=v\left( t \right)+\frac{\delta t}{2}\left[ a\left( t+\delta t \right)+a\left( t \right) \right]+O\left( \delta {{t}^{3}} \right)$

速度Verlet算法具有以下优点：可以同时进行坐标r(t＋δt)和速度v(t＋δt)计算，并且精度也相同；适合(r_N，v_N)轨迹作用的粒子体系.

2.3 蛙跳算法(Leap-frog algorithm)

Hockney提出的蛙跳算法^[12]，是由Verlet算法演变而得到的算法. 蛙跳算法相比Verlet算法有2个优势：①它带有显速度项；②它无需计算两个较大的量2r(t)和r(t-δt)的差，减小计算量的同时也减小误差.但它也有一个明显的缺点：位置坐标r(t+δt)和速度v(t+$\frac{1}{2}$δt)不能同步，即系统中各个粒子的位置确定下来之后，系统的势能才随之确定，但却不能计算出这个时刻系统的动能对总能量的贡献.

2.4 Beeman算法

Beeman算法^[13]是由Beeman提出的一种与Verlet算法相类似的算法，是另一种运用较为广泛的算法：

$\begin{align} & {{r}_{i}}\left( t+\delta t \right)={{r}_{i}}\left( t \right)+{{v}_{i}}\left( t \right)\delta t+\frac{1}{6}\left[ 4{{a}_{i}}\left( t \right)-{{a}_{i}}\left( t-\delta t \right) \right]\delta {{t}^{2}} \\ & {{v}_{i}}\left( t+\delta t \right)={{v}_{i}}\left( t \right)+\frac{1}{6}\left[ 2{{a}_{i}}\left( t+\delta t \right)+5{{a}_{i}}\left( t \right)-{{a}_{i}}\left( t-\delta t \right) \right]\delta t \\ \end{align}$

其优点是在准确性一致时可以引用较长积分步长；缺点是储存量较大.

2.5 Rahman算法

Rahman算法^[14]是Rahman首次采用的一种积分方法，它属于预测—校正算法的变通形式，尽管Rahman算法求解运动方程能得到较为准确的解，但缺点是计算量大，耗时长，现在使用很少.

3 边界条件

即使是现代的巨型计算机，MD方法还是只能用于模拟几百到几千粒子数的体系，这就值得我们思考：如何用少量粒子体系来模拟宏观体系，周期性边界条件的引入很好地解决了这个问题^[15]. 周期性边界条件为各个方向上原胞的影像提供周期性的排列，使得我们可以模拟少量粒子的体系来研究宏观物质的宏观性质. 在体系中通常采用“最小影像”约定处理粒子与粒子之间的相互作用，即一个粒子与其他无穷重复的MD模拟基本计算原胞中的另外N-1个(假设基本原胞中有N个粒子)粒子中的每个粒子或其最邻近影像粒子发生相互作用. 但并非所有体系都适用于周期性边界条件，有必须用到周期性边界条件的体系，当然也有必须用到非周期性边界条件的体系，比如自带边界条件的液滴和团簇；也有一些可能需要用到非周期性边界条件的体系，如非均匀体系或处在非平衡的系统.

4 势函数

MD模拟过程中的粒子因为受力而不停的运动，原子受到的力与其所在的空间位置有着直接的关系，用于表达两者的密切关系具体数学公式便是势函数.体系中任何原子的运动都受到势函数的约束，可以说势函数引导着体系微观结构的构成，从而决定模拟对象的宏观性质. 所以说在MD模拟中，势函数的选择直接决定着模拟结果的正确与否. 目前,分子动力学模拟技术还处于发展时期,存在着势函数数量有限、模拟规模较小、耗时长、模拟模型较简单等问题^[16]. 原子之间的作用势一直是MD模拟的重点难点，但其缓慢的发展在一定程度上制约了MD方法在实际研究中的应用.原子间的相互作用势从最初的对势发展到如今的多体势，对势认为原子间相互作用只存在两两原子之间，与其他原子的位置无关，也不考虑局域内电子的作用效应；多体势发展为其势能不局限于两两原子核之间，原子周边的电子所带有的能量也是体系势能的一部分. 而对于实际多原子体系，一个原子的位置变化将会影响其他原子间的有效相互作用，因此，多体势更适合多原子体系.

5 分子动力学模拟的系综

系综是采用统计力学描述系统的统计规律性时引入的一个最基本的概念，它是指在一定宏观条件下，大量性质和结构完全相同的系统的集合. 分子动力学引入系综这个概念具有非常重要的意义，因为所有的MD模拟必须在有系综的前提下才能开始.常用的系综有以下3种：

1） 正则系综，全称为宏观正则系综，简写为NVT.即保持体系的粒子数N，体积V和温度T不变.通常采用热浴法和速度标定法来保持温度的恒定，其中热浴法是NVT系综中最常用的方法，简单来说就是将系统与外界的热源处于热力学平衡状态，而速度标定法则是直接通过标定速度来实现温度的恒定.

2） 微正则系综，简写为NVE，顾名思义就是保持体系中的粒子数N，体积V和总能量E的不变.NVE系综下的体系被看作是一个孤立的平衡态，与外界没有粒子的交换，也没有能量的交换，体系中的势能和动能随着体系的运动而发生改变，但是它们的总和不会发生改变.通常体系的初始模型的能量是不能被确定的，为了把能量调整至合适的值，必须先进行一定时间的模拟，根据能量的变化进行调节从而确定体系的固定能量值.但在这个过程中，体系可能需要打破已达到的平衡状态，所以必须给予足够的时间满足体系再次建立平衡.

3） 等温等压系综，简写为NPT，即表示体系具有确定的粒子数N，压强P和温度T.温度的恒定是通过调节系统的速度来实现，然而由于体积V和压力P是共轭量，所以调节压力的恒定更加复杂，目前可采用标度体系的体积和外部热浴法来实现压力的恒定.

6 热力学性质的计算

在模拟得出体系粒子在相空间的运动轨迹后，运用物理统计学原理即可得到体系的结构特点和宏观特性，就冶金熔体的结构分析而言，目前常用的分析方法有：径向分布函数、均方根位移和配位数.

6.1 径向分布函数(RDF)

RDF表征的是结构的无序化程度，其表达式为：

$g\left( r \right)=V/N\cdot n\left( r \right)/4\pi {{r}^{2}}\Delta r$

其中，n(r)表示距离原点r到r+Δr之间的平均原子数，模拟体系的体积为V.它所描述的是一个原子周围其他原子的分布情况，反映的是在一个粒子周围(距离r)发现另一个粒子的概率.

6.2 均方根位移

在MD模拟计算中，体系粒子从计算开始就一直在运动，位置也一直在发生变化.这里把粒子位移平方的平均值定义为均方位移，简写为MSD，其表达公式为：

$MSD=\left\langle {{\Delta }^{2}}r\left( t \right) \right\rangle =\frac{1}{N}\sum\limits_{i=1}^{N}{{}}\left\langle {{\left[ {{r}_{i}}\left( t \right)-{{r}_{i}}\left( 0 \right) \right]}^{2}} \right\rangle $

r_i(0)和r_i(t)分别表示原子i在零时刻和t时刻的位移.根据爱因斯坦的扩散定律，当体系处于液态时，其原子的扩散行为可用均方位移随时间的变化来表示，关系式为$\underset{r\to \infty }{\mathop{\lim }}\,$MSD=c+6Dt. 其中D是扩散系数，c是常数.由公式可看出均方位移的大小与原子的活跃度有关，均方位移越大则表示体系中原子越活跃，体系中的扩散系数也越大.

6.3 配位数

配位数就是原子的第一近邻原子的个数.配位数的多少表示该原子周围的原子分布的密度大小,它经常在结构分析中用作为一种辅助手段^[17].

7 冶金熔渣结构的分子动力学模拟

炉渣又称为熔渣，是炼钢过程中产生的，在实际生产中常说：炼钢就是炼渣，所以炼好渣是实现炼钢生产优质、高产、低消耗的重要保证，为此有必要更加深入地了解炉渣的结构性质. 很多研究者^[18-21]用常规实验方法对一些冶金炉渣的结构模型进行热力学性质的研究，并且也取得很好的进展，但由于熔渣所需的高温环境，难以通过常规实验来更加深入地了解炼钢过程中熔融炉渣的性质，MD模拟的引入恰好可以很好地解决这个问题，而且还能通过计算机图像功能展示熔渣的微观结构. 目前在冶金炉渣领域的应用还较少，早些年W.Smith等^[22]就对Na₂SiO₅进行分子动力学模拟，尽管结果分析表明非桥氧的相对分布与核磁共振数据有所差异，但这也为未来的模型计算提供一个非常严格的标准；Sasaki Y等^[23-24]为了确定F离子对熔融碱金属硅酸盐(Na₂O-NaF-SiO₂、Na₂O-NaF-Al₂O₃-SiO₂)体系影响，对其进行分子动力学模拟，证实F取代O增强了体系(Na₂O-NaF-SiO₂)的聚合度,而在体系(Na₂O-NaF-Al₂O₃-SiO₂)中随着NaF的增加铝氧四面体增加，证实NaF可以当做铝四面体连接的网格修饰子，但对Si四面体没用；张毅刚等^[25]运用MD方法研究Mg₂SiO₄在高温高压下条件下的黏滞度，在压力为135 GPa和温度达到 5 595 K时，熔体的黏滞度与温度成反比；邓莉等^[26]对高温高压下的MgSiO₃熔体做了相应的研究，他们指出随着压力的不断增加，原子间的平均配位数随之增加，熔体结构变得更加致密，而Si-O多面体结构逐渐过渡为相对稳定的Si-O八面体结构；吴永全等^[27-31]利用MD方法对硅酸盐炉渣体系(CaO-Al₂O₃、CaO-SiO₂、Al₂O₃-SiO₂和CaO-Al₂O₃-SiO₂)做了一系列的研究，并且与高温拉曼检测结果做了对比分析. 研究表明，常压下Si、Al在炉渣中主要以SiO₄^-和AlO₄^-四面体的结构形式存在，而且SiO₄^- 四面体稳定性比AlO₄^-四面体要好. Si、Al在硅酸盐熔体中主要起网格形成子的作用，而Ca是典型的网格修饰子，主要起破坏网格的作用，氧除了存在于桥氧、非桥氧中，还发现少量不稳定的自由氧，并且还得到3类氧之间的平衡关系式. 除此之外还计算得到了SiO₄^-四面体结构中Si-O键长和键角等有关数据；肖建平^[32]对不同组成成分的钒渣体系(FeO-V₂O₃、FeO-V₂O₃-SiO₂和FeO-V₂O₃-TiO₂)进行研究，结果表明随着温度的降低FeO-V₂O₃体系中发现了长程有序结构.对体系FeO-V₂O₃、FeO-V₂O₃-TiO₂和FeO-V₂O₃-SiO₂中的键角和配位数的结合分析可得Fe-O呈四面体配位，V-O为八面体配位.另外还得出往FeO-V₂O₃渣系中加入SiO₂会使得粒子的扩散系数变大；张燮^[33]运用第一性原理、Mulliken电荷分析以及曲线拟合的方法对FeO-SiO₂-V₂O₃三元渣建立一套新的势参数，对体系在不同组成和温度下进行MD模拟，并且通过统计力学方法建立渣中各离子和熔渣黏度与结构之间的关系，从而实现对渣中离子系数与剪切黏度的计算；张国莹等^[34]在固定成分Al₂O₃/SiO₂=3/2的情况下，利用MD模拟不同比例含量的Na₂O对Na₂O-Al₂O₃-SiO₂体系结构的影响，结果表明当Na⁺含量较高时，其本身可以当作起解聚网络作用的网络修饰子；王明国^[35]运用MD模拟Fe₂O₃含量的变化对Ca:Al:Si比例为1:1:3的CaO-Al₂O₃-SiO₂-Fe₂O₃体系结构的影响，研究表明体系中的Fe³⁺难以形成[FeO₄]结构存在于网络结构中，因为在与O₂-的结合能力上Si⁴⁺和Al³⁺要比Fe³⁺更强.另外，在体系中还发现三配位氧的存在，其产生主要是因为体系中的阳离子对O₂-的争夺；Ting WU等^[36]还利用分子动力学模拟方法研究CaO-SiO₂和CaO-Al₂O₃体系微观结构和运输性质，并分析两者之间的关系.研究表明CaO在体系中起解聚作用，随着CaO含量的增加，熔体的结构被解聚，使得熔体的黏度降低和电导率增强.

从上述研究证实MD模拟在冶金炉渣结构的研究是很有意义的，通过分子动力学计算机模拟，越来越多的微观结构信息被挖掘，同时对微观结构性质了解也有助于我们对宏观性质的理解.尽管目前MD模拟在冶金炉渣方面的应用比较少，但随着钢铁冶金技术和计算机模拟的发展，两者的交集必然会越来越多.

8 冶金炉渣的分子动力学模拟发展方向

分子动力学模拟在冶金炉渣结构中应用已经取得不错的进展，研究从微观角度描述熔渣的结构演变，从而构建沟通炉渣微观和宏观性质的桥梁.在分子动力学模拟过程中，最重要的因素是体系的初始构象以及势函数的选取，初始构象由前人对熔渣体系的理解以及大量的模拟实验给出，而势函数是整个模拟的关键，它的选取直接决定着模拟结果的准确性，不同的势函数也有不同的适用范围.在势参数上，目前大多势参数还是采用经验参数，然而非经验参数的准确性比经验参数有优势，但它在拟合计算时比较麻烦.所以找到一个在冶金熔渣上模拟精确又不复杂的势函数将会是一个研究的热点.

就目前的炉渣结构模拟而言，体系还是相对比较简单.但实际炉渣的组成和结构远比目前研究的要复杂，很多稀少但活跃的元素并没有被列为研究的对象，但它们的存在对渣的性能有极大的影响.另一方面，目前的分子动力学模拟技术只能模拟几百到几千个原子的体系，这对那些多组成成分的复杂体系的计算也将是个挑战.故分子动力学计算机模拟技术的突破和多组成成分的复杂体系的模拟研究也将是从事炉渣模拟研究工作者的一个重要的研究方向.