1. 引言
20世纪70年代，Aisenberg、Pencen等人首次利用离子束沉积法(IBD)在室温下制得了一种硬质膜，称其为类金刚石(DLC)薄膜[1] [2] 。到70年代末期，E. G. Spencer和Weissmantel等提出了离子束增强沉积法(IBED)制备DLC膜，并率先对其进行了研究[3] 。由于制备方法简单，同时具有金刚石和石墨的优点,以及在光学，电学，机械制造，医学，惯性约束聚变等领域拥有广阔的应用前景，这种类金刚石薄膜在科研界迅速引起了极大的关注。
随着不断的研究人们已经逐渐认识到这种类金刚石薄膜是一种非晶薄膜，拥有网状的结构。并发现非晶类金刚石薄膜同时含有sp3、sp2、sp1这三种杂化键，它的性能随sp3、sp2、sp1的含量不同而发生变化[4] 。在信息时代，为了更好的研究性能和结构，仅仅依靠传统的实验方法显然已无法满足需要，而分子动力学模拟在材料研究中所表现出的省时，省料，以及节约劳动力的巨大优势，使这种方法得到了迅速发展。所谓分子动力学模拟是指以玻恩-奥本海默近似[5] 为基本前提，选取合适的系综和原子相互作用势建立模型，并用牛顿力学方程对体系中的粒子进行求解，直到体系性质趋于稳定。可以说实验与模拟相结合的研究方法越来越受研究者们的青睐，它实现了人们对材料有目的、有规律、高效率的进行研究的梦想。
2. 类金刚石薄膜的分子动力学研究
2.1. 原子作用势函数
上世纪80年代，美国的D. Beeman和英国的J. Robertson等人对无序的非晶碳模型进行了研究，发现体系中的sp2结构所得的数据与实验结果差距较大，这是因为在模拟中原子的相互作用势对π键的考虑不足造成的[6] -[8] 。紧接着，IBM研究部门的J. Tersoff引入了一种经验原子势函数，研究发现，这种势函数可以非常精确的描述碳的结构、性能等[9] [10] 。但Tersoff势在描述π-π键交互能上依然无能为力。1990年，美国的Brenner基于Tersoff势引入了经验的多体势函数，由于在计算上具有短程和快速的优点，这使得进行大规模的分子动力学模拟得以实现[11] ，但其仍不能对石墨的结构进行完整的模拟。后来，克里特岛大学的Pantelis C Kelires用修正后的势函数[12] ，对非晶碳模型的结构，力学性能和弹性性能进行了研究。发现模型的弹性软化和局部应力对网状非晶碳的形成起了促进作用，但是对sp3到sp2的转化机制和电子的行为并没有给出很好的解释[13] [14] 。
为了更好的考虑电子行为对类金刚石薄膜性能的分子动力学模拟的影响，爱荷华州立大学的Wang等人采用紧束缚的分子动力学方法计算了类金刚石薄膜电子结构，但遗憾的发现高密度的sp3含量与试验结果存在很大的偏差[15] [16] 。德国的Th. Frauenheim等人则用半经验势密度泛函近似的分子动力学对淬火的非晶碳模型和加氢的非晶碳模型的结构，物理性能和电子性能进行了模拟研究，并利用二体交换势和局域态密度近似将多电子相互交换能等效为一个电子的势能。发现这可以对π键的影响给出很好的解释，并指出高密度，高四面体键的模型的电子带性能与sp3模型中嵌入的π键的尺寸和重叠分布有一定的关系，p杂化轨道的重叠促进了强π键的形成[17] -[20] 。后续还有许多人做了相关方面的研究，他们从新的模型结构，大的计算基组，以及第一性原理分子动力学等方面入手，寻求更高的精确度，在前人的工作基础上取得了一定的进展[21] -[29] 。
我们可以发现，采用Tersoff势和Brenner势可满足对具有大量粒子、且电子行为与π-π键交互作用对其性质影响不大的系统的分子动力学研究的需要。但当π键与电子行为等因素对薄膜性质有较为重要的影响时，上述两种作用势将不能满足对类金刚石薄膜分子动力学研究的要求，在这种情况下，为了更好的逼近真实情况，并得到更好的计算精度，更加准确的原子间作用势是必须考虑解决的最重要的影响因素之一。例如，改进后的多体Tersoff势函数，甚至是采用精度更高的第一性原理法等都是目前研究者们所做的一些有益的尝试[17] [23] [24] 。但另一方面，我们也应该认识到，随着计算精度的提高，分子动力学模拟对计算的软硬件的要求也越来越高，而且，也会大幅增加计算量和计算时间，因此，在实际的模拟工作中，研究者必须根据研究的实际需要，采用合适的模型和原子作用势来进行计算，以获得计算精度与计算工作量之间的最佳效能比。
2.2. 国内外对薄膜生长机理的研究
到目前为止，针对类金刚石薄膜的生长机理已有很多理论与实验研究工作，但还未取得一致的认识，比较获得接受的生长机理是认为，由浅层碳离子注入所引起的内部亚表面原子的生长[30] 。例如，德国科学家H. U. Jager和A. Yu. Belov以Brenner势模拟原子间的相互作用，用分子动力学模拟研究了ta-C薄膜的形成，并进一步研究了薄膜在退火处理后网状非晶碳中sp2团簇结构的形成与演化[30] -[32] 。他们发现，基体的温度和沉积粒子的入射能量对薄膜的生长有显著的影响，具体表现为：在低温下薄膜中sp3键含量随能量的增加而增加，但当沉积离子束的能量低于20 eV时，将不会产生具有类金刚石性能的薄膜，这与Y. Lifshitz等提出的亚表层注入模型的情况相一致[33] -[35] 。当基体的温度超过临界温度Tc(约为100℃)时，ta-C将急剧转变为石墨相。在对薄膜进行后续的退火处理研究时，通过比较退火薄膜与沉积型薄膜的结构发现，当薄膜的密度小于3 g/cm3时，随着大量的sp2键的形成，ta-C的结构渐趋不稳定，这与早期的相关报道也是符合的[30] -[32] 。但是，由于研究中采用的原子间作用势是Brenner势函数，正如我们上文中所讨论的，由于该势无法很好的模拟π键的影响，因此，应用该势很难对富含sp2键的薄膜性质作出很好的预测，这方面还有待于后续研究，例如构建更精确的原子间作用势等，以及与实验研究的结合，才能对富含sp2键的类金刚石薄膜的结构与性质有一个更好的认识。
在类金刚石薄膜生长机理的分子动力学研究方面国内学者也开展了大量的工作，例如，清华大学马天宝等人以半经验的反应式键序势函数(REBO) [36] ，用分子动力学从原子尺度上研究了薄膜的生长[37] 。他们的研究指出，原子在薄膜表面的行为可分为表面冷冻，注入，迁移和反弹等四种，并给出了四种行为与能量变化的关系，如图1所示。研究结果显示，在低能量态下，原子主要存在迁移和冷冻行为，当能量超过约为30 eV的阀值时，原子的注入行为随着入射能量的升高显著增加，而冷冻行为显著减少。在高能量态下，原子的行为主要为注入和迁移。原子注入易于形成sp3杂化占优的致密薄膜，而原子迁移却导致薄膜的热弛豫，使亚稳态sp3结构向更稳定的sp2结构转化。入射原子注入和迁移数目及剧烈程度的对比决定了DLC薄膜生长的模式和最终结构，即类金刚石薄膜的最终结构是由于原子注入效应和
Figure 1. Relationship between surface atomic behavior and the incident energy [37]
图1. 表面原子行为与入射能量的关系[37]
弛豫效应互相竞争的结果[37] 。马天宝等人的研究解决了Lifshitz的亚表层注入模型在理论上存在的局限性，具体的阐述了薄膜表面的原子弛豫现象，对研究入射原子的表面行为具有重要贡献，并且对sp3占优的薄膜的制备方法提出了有价值的理论建议。
2.3. 类金刚石薄膜性能的研究
国外关于类金刚石薄膜的摩擦性能的研究早在上世纪80年代就已经开展[38] ，但在国内起步的相对较晚。王丽莉等人用基于密度泛函理论的第一原理分子动力学方法模拟研究了类金刚石薄膜(110)面的结合强度和摩擦性能[39] 。研究发现，表面排斥相互作用发生在0.25 nm以内，表面的结合强度为6 Mpa。而通过<001>和<110>两个方向上的摩擦性能的对比，发现<110>方向上的摩擦力一直高于<001>方向上的摩擦力，且<110>方向上的滑动磨损率更高，更易于发生摩擦表层的偏转。该项研究分析了摩擦力各向异性的成因，从理论上对类金刚石薄膜抛光率的各向异性做出了定性的解释[39] 。
类金刚石薄膜与基体的结合力是影响类金刚石薄膜应用的重要性能的之一，理论研究发现，其可以通过多叠层结构来增强[40] -[42] 。E. B. Halac等人利用分子动力学方法，通过改变沉积能量的方法得到了多叠层结构的非晶碳薄膜模型，并对多叠层结构对膜基结合力的影响进行了细致的研究[43] 。他们的研究发现，粒子入射能量的增大有益于提高薄膜的抗压强度、密度、sp3含量等性质，而多叠层结构则可以降低薄膜的内应力，进而增强膜基结合力和薄膜的热稳定性。并且，通过对薄膜进一步退火研究发现，多叠层样品只有最外面的浅层区域发生了结合键的转变和抗压应力降低的现象，而大部分内部区域的sp3/sp2比率几乎不变，这说明多叠层结构的热稳定性好，且证实了非晶碳薄膜的多叠层结构很好的保护了退火后样品的硬度等性能。该项研究表明，可以通过改变能量来获得不同参数的多叠层结构的薄膜，多叠层结构利于提高膜基的结合力，但随层数的增加，薄膜机械和电子性能的演变则还需要更加深入的研究。
此外，通过分子动力学的方法还可以模拟不同的工艺条件对类金刚石薄膜表面粗糙度的影响，例如，Minwoong Joe等人通过分子动力学研究指出，类金刚石薄膜表面粗糙度与入射原子团剂量及其入射角度均有很大的依赖关系，计算结果如图2所示。在入射角度为0˚和30˚时，薄膜的粗糙度几乎不随剂量的变化而变化；当入射角度变为60˚和70˚时，薄膜的粗糙度随入射剂量的增加显著增大。他们经过综合分析后发现，随着入射角度的增大，入射粒子的冲击会诱导薄膜表面原子向上传输以形成晶核结构，进而对
Figure 2. Surface roughness evolution as a function of ion dose. Three stages of roughness evolution are shown: the initial roughening stage (I), the seed-forming stage (II), and the shadowing-driven rapid roughening stage (III) [44]
图2. 表面粗糙度随离子剂量的演化图。粗糙度的演化分为三个阶段展示：最初粗糙度阶段(Ⅰ)，晶核形成阶段(Ⅱ)，遮蔽引起的快速粗化阶段(III)[44]
后续入射原子团产生遮蔽效应，最终导致薄膜表面粗糙度的急剧增大[44] 。尽管他们的研究结果对如何控制薄膜表面粗糙度提出了可能的途径，但通过对实验研究的调研可以发现，粒子的入射能量也是造成薄膜表面粗糙度的重要影响因素之一，而在Joe等的研究中并没有考虑入射能量的影响。
经过二、三十年的发展，除了上文中涉及的针对薄膜生长及其物理性质的模拟外，利用分子动力学方法对类金刚石薄膜的研究工作仍然方兴未艾，目前已拓展到诸多领域，例如，关于新型沉积源的模拟研究[45] 、以及和纳米学，细菌学，电学等的结合等[46] -[49] 。这充分表明，针对类金刚石薄膜的分子动力学研究已发展为一种面向应用的多学科综合性研究。
3. 结束语
本文较为系统的综述了针对类金刚石薄膜的分子动力学研究的进展。以往的研究表明，原子间相互作用势是影响分子动力学模拟计算结果准确性的最重要的因素之一。Tersoff势和Brenner势完全可以满足粒子数量巨大、但电子行为与π-π键交互作用对其性质影响不大且对计算精度要求不高的类金刚石薄膜的分子动力学研究。但若需得到高精度的计算结果，则需考虑采用大基组或电子行为与π键影响的多体势的分子动力学方法，或者采用第一性原理的分子动力学方法也是一种有益的尝试[23] [24] [27] 。此外，尽管针对类金刚石薄膜性能模拟与优化的分子动力学研究已有很多，但通过对当前研究结果的综合分析后不难发现，在对类金刚石薄膜的生长机理，以及对其膜基结合力、表面粗糙度、热稳定性等物理性质的认识和优化方面还存在很多尚未解决的问题，这都有待于我们进一步的深入研究。
致 谢
感谢国家重大专项和北京工业大学先进技术基金的资助。
NOTES
*通讯作者。