多尺度模拟方法在材料位移损伤 效应研究中的应用

贺朝会, 唐 杜, 李 奎

(西安交通大学 核科学与技术学院, 西安 710049)

辐射在材料中产生位移损伤是一个从微观到宏观的过程，一般经历离位级联、热峰阶段、缺陷初步演化和缺陷迁移4个阶段[1]，最终导致材料宏观性能的改变。实验中可观察的是宏观性能的变化过程，微观变化过程只有通过数值模拟才能揭示。位移损伤效应的模拟难度很大，因为它涉及一个多尺度的变化过程，材料尺寸从亚纳米到米量级，时间从亚皮秒到年量级，如图1所示，这需要用到多种不同的方法和软件进行多尺度的模拟研究[2]，如图2所示。

图1 材料性能多尺度变化示意图Fig.1 Multi-scale sketch map of material performance change

图2 多尺度数值模拟方法示意图Fig.2 Multi-scale sketch map of numerical simulation methods

1 第一性原理[3-9]

第一性原理方法(first-principle method)也被称为从头算方法(ab initio methods)，是指基于量子力学理论，完全由理论推导而得，不使用除基本物理常数和原子量以外的实验数据以及经验或者半经验参数的求解薛定谔方程的计算方法。基本思想是：将多原子构成的体系理解为电子和原子核组成的多粒子系统，根据原子核和电子互相作用的原理及其基本运动规律，运用量子力学基本原理最大限度地对问题进行“非经验性”处理。

量子力学描述n粒子体系的波函数包含3n个坐标，相应的薛定谔方程是含3n个变量的偏微分方程，而针对绝大部分问题涉及的体系，这就会是一个非常复杂甚至无法解决的问题。

密度泛函理论(density-functional theory，DFT)采用电子密度分布n(r)作为基本变量，而非波函数，研究原子、分子、固体等多电子体系的基态性质。DFT建立在Hohenberg和Kohn提出的非均匀电子气理论基础之上，可以归结为两个基本定理。定理1：不计自旋的全同费米子系统非简并基态的所有性质都是粒子密度函数的唯一泛函。该定理保证了粒子密度作为体系基本物理量的合法性，同时也是DFT名称的由来。定理2：DFT的变分法中，对于一个给定的外势，真实电子密度使能量泛函取得最小值。由于DFT中融入了统计的思想，不必考虑每个电子的行为，只需计算总的电子密度，计算量大为减少，从而大大简化了应用量子力学计算材料性质所涉及的复杂数学问题。DFT把多电子问题成功转化为单电子问题，近似求解薛定谔方程，可用于描述原子、分子等微观粒子体系的量子行为，分析其电子结构，计算总能量。已广泛应用在材料设计、合成、模拟计算和评价等诸多方面，成为计算凝聚态物理、计算材料科学和量子化学的重要基础和核心技术。

广义的第一性原理包括两大类，以Hartree-Fork自洽场计算为基础的ab initio从头算和DFT计算。也有人主张，ab initio专指从头算，而第一性原理和所谓量子化学计算特指DFT计算。

第一性原理主要应用于1)材料物性计算，包括晶格常数、弹性模量、相变、电介常数、热膨胀系数、光学特性等；2)电子结构计算、力学行为预测；3)化学成键机理描述；4)半导体材料能带计算，包括缺陷形成能与结合能、晶格原子的离位阈能、杂质或缺陷的迁移能；5)分子、团簇、晶体等的结构优化及新构相预测；6)纳米器件设计；7)极端条件下的材料行为模拟等。例如，CdSxTe1-x三元半导体材料的基态物性计算，硅材料中的间隙缺陷结构及扩散路径，硅材料中不同电荷态空位和间隙原子的形成能，半导体材料AlSb点缺陷迁移率退化等。

2 分子动力学[1-14]

分子动力学(molecular dynamics，MD)是一种确定性方法，按照体系内部的内禀动力学规律来确定位形的转变，跟踪系统中每个粒子的个体运动。根据统计物理规律，给出分子的坐标、速度等微观量与温度、压力、比热容、弹性模量等宏观可观测量之间的关系。

对于多体系问题的严格求解，需要建立并求解体系的薛定谔方程，根据波恩-奥本海默近似，原子核的运动可以用经典动力学方法处理。把薛定谔方程的求解简化为牛顿运动方程的计算。具体步骤包括：1)建立一个由n个粒子(分子)组成的模型体系；2)解n个粒子(分子)组成的模型体系的牛顿运动方程直至平衡；3)平衡后，进行材料性能的计算，对模拟结果进行分析。

MD方法只考虑多体系中原子核的运动，而不考虑电子的运动，忽略量子效应。适用于经典近似，在很多材料体系应用中都较精确。不适用于涉及电荷重新分布的化学反应、键的形成与断裂、解离、极化以及金属离子的化学键等领域。不适用的领域可用量子力学方法替代，也有结合经典分子动力学和量子力学的新方法，叫做第一性原理分子动力学。

进行MD运算有7个必备步骤：1)建立计算模型；2)设定计算模型的初始坐标和初始速度；3)选定合适的时间步长；4)选取合适的原子间相互作用势函数，便于进行力的计算；5)选择合适的算法、边界条件和外界条件；6)计算；7)对计算数据进行统计处理。MD运算主要用于研究材料微观性能。如单晶硅辐照级联碰撞过程分子动力学模拟。图3为初级反冲原子(primary knock-out atom, PKA)产生的缺陷及其随时间演化的过程。图4为2 keV 硅PKA产生的点缺陷的数目随时间的变化情况。主要经历3个阶段：初始阶段——缺陷数目急剧上升；热峰阶段——缺陷数目达到顶峰并开始下降；趋衡阶段——缺陷数目和构型趋于稳定。

3 动力学蒙特卡罗[14-20]

MD方法可以模拟材料中产生位移损伤的离位级联、热峰阶段和缺陷初步演化的前3个阶段，时间尺度到纳秒量级、尺寸尺度达亚微米量级。而对于第4阶段的缺陷迁移的模拟，MD方法则无能为力，只有依靠动力学蒙特卡罗(dynamics Monte Carlo，KMC)方法。KMC方法是将普通的蒙特卡罗方法中的事件赋予时间的尺度，模拟系统随时间演化的过程。可以模拟界面、可动粒子、扩展缺陷、缺陷团簇和复杂团簇的迁移、扩散、复合、发射粒子、捕获和解离等，时间尺度达105s。模拟中所有事件的发生频率遵从Arrhenius方程：

vij=v0·exp(-Eij/kBT)

(1)

其中,v0为指前因子；Eij为某粒子从状态i转变为状态j所需越过的势垒；kB为玻耳兹曼常数；T为温度。图5为KMC模拟241.9 keV硅PKA产生的不同缺陷随时间的演化。缺陷退火分3个阶段：Ⅰ.1×10-11s ≤t< 2×10-3s：点缺陷成团；Ⅱ.2×10-3s ≤t< 2×102s：缺陷团内间隙原子和空位发生复合反应；Ⅲ.t≥ 2×102s：小缺陷团缓慢发射点缺陷。

4 待研究问题

虽然已经开展了一些研究工作，但仍然存在许多亟待研究的问题。

1) 基本参数研究

在应用LAMMPS进行MD模拟计算时发现一些材料的势函数还存在问题，如GaAs，SiGe等；应用KMC进行缺陷演变规律研究，模拟的缺陷种类较少，一些类型的缺陷团的形成和解离反应还无法模拟，需添加各种缺陷反应类型和反应参数。

2) 不同辐射源产生的缺陷演化规律异同性研究

图3 PKA产生的缺陷及其随时间演化的过程Fig.3 Defects produced by primary knock-out atom and their evolution with time

图4 2 keV Si PKA产生点缺陷的数目随时间的变化情况Fig.4 The evolution of the number of point defects produced by 2 keV Si PKA with time

(a) The evolution of interstitials and vacacies with time

(b) The evolution of two,three and four vacacies with time

(c) The evolution of impurity related defects with time 图5 241.9 keV硅PKA产生的不同缺陷随时间的演化Fig.5 The evolution of the different defects produced by 241.9 keV Si PKA with time

针对不同辐射源产生的位移损伤，已有初步的等效性研究结果，主要是研究材料的宏观性能的改变，其深层次的原因有待揭示。基于缺陷的多尺度模拟方法，可以研究中子、质子、重离子和电子产生的缺陷及其演化规律的异同和对Si器件性能影响的差异；从缺陷的产生和演变的角度分析不同辐射源之间的等效性，揭示其内在机理；研究均匀分布缺陷和非均匀分布缺陷以及不同浓度的缺陷演化规律的异同。

3) 单粒子位移损伤研究

已经开始了Si器件的单粒子位移损伤研究，可以进一步开展二元半导体材料，如SiC、GaAs、GaN、SiGe等的单粒子位移损伤研究、光电材料和器件的单粒子位移损伤研究，以及半导体辐射探测器的位移损伤漏电流研究。

4) 位移损伤缺陷长时间演化的温度效应研究

由于不同温度下缺陷会有不同的演化行为，例如V+O在室温下形成VO后是稳定的，不会进一步形成V2O和VO2缺陷，但当温度高于300 ℃时，可能形成V2O和VO2两种缺陷。因此，可以利用MD模拟PKA在不同温度下的缺陷初态分布；并将不同温度下的缺陷初态分布引入KMC模拟中，进一步研究不同温度下各种缺陷的长时间演化行为。

5) 软件研究

要开展多尺度的模拟研究，不同软件之间的衔接是个问题，目前仅实现了MD结果到KMC输入的转换，如图6所示，而其他软件之间的衔接难度和工作量都很大，且MD结果到KMC输入的转换目前是通过独立编程的软件实现，还没有实现自动化。此外，软件的国产化更是一个长期、艰巨的任务，可以从编制用于二元半导体材料如GaAs、GaN、SiC、SiGe等的KMC软件入手，首先实现KMC软件的国产化。

实际要研究的问题还有很多，需要进一步细化分析。

(a) The introduction of the number and spatial distributionof defects from MD to KMC

(b) The evolution of the number of defects with time 图6 MD与KMC的衔接Fig.6 The connection between MD and KMC

5 结论

应用第一性原理、分子动力学和动力学蒙特卡罗等方法，可以跨尺度，尺寸从亚纳米到米量级，时间从亚皮秒到年量级，模拟材料中的位移损伤效应，揭示辐照效应机理。多尺度数值模拟方法为材料辐照效应研究提供了有力的工具，但仍然存在许多问题亟待研究。

[1]NORDLUND K, DJURABEKOVA F. Multiscale modelling of irradiation in nanostructures[J]. Journal of Computational Electronics, 2014, 13 (1): 122-141.

[2]DE LA RUBIA T D, CATURLA M J, ALONSO E A, et al. The primary damage state and its evolution over multiple length and time scales: Recent atomic-scale computer simulation studies[J]. Radiation Effects and Defects in Solids, 1999, 148 (1): 95-126.

[3]KRESSE G, HAFNER J. First-principles study of the adsorption of atomic H on Ni (111), (100) and (110)[J]. Surf Sci, 2000, 459(3): 287-302．

[4]JI Y F, WANG B, LUO Y. First principles study of O2adsorption on reduced rutile TiO2-(110) surface under UV illumination and its role on CO oxidation[J]. J Phys Chem C, 2013, 117(2): 956-961.

[5]SEBBARI K, ROQUES J, SIMONI E, et al. First-principles molecular dynamics simulations of uranyl ion interaction at the water/rutile TiO2(110) interface[J]. Surf Sci, 2012, 606 (15/16): 1 135-1 141.

[6]BENDAVID L I, CARTER E A. First-principles predictions of the structure, stability and photocatalytic potential of Cu2O surfaces[J]. J Phys Chem B, 2013, 117 (49): 15 750-15 760.

[7]ISLAM M M, DIAWARA B, MAURICE V, et al. First principles investigation on the stabilization mechanisms of the polar copper terminated Cu2O(111) surface[J]. Surf Sci, 2009, 603 (13): 2 087-2 095.

[8]LI K, ZHAO Y L, DENG J, et al. Adsorption of radioiodine on Cu2O surfaces: A first-principles density functional study[J]. Acta Phys Chim Sin, 2016, 32 (9): 2 264-2 270.

[9] 李奎. 碘分子、碘阴离子在Cu2O和CuFeS2表面吸附的机理研究[D]. 西安: 西安交通大学, 2016. (LI Kui. Study on adsorption mechanisms of iodine anions on Cu2O and CuFeS2surfaces[D]. Xi’an: Xi’an Jiaotong University, 2016.)

[10] PLIMPTON S. Fast parallel algorithms for short-range molecular-dynamics[J]. Journal of Computational Physics, 1995, 117 (1): 1-19.

[11] BORODIN V A. Molecular dynamics simulation of annealing of post-ballistic cascade remnants in silicon[J]. Nucl Instrum Methods Phys Res B, 2012, 282: 33-37.

[12] CATURLA M J, DE LA RUBIA T D, GILMER G H. Disordering and defect production in silicon by keV ion irradiation studied by molecular dynamics[J]. Nucl Instrum Methods Phys Res B, 1995, 106 (1/4): 1-8.

[13] SANTOS I, MARQUÉS LA, PELAZ L, et al. Temperature effect on damage generation mechanisms during ion implantation in Si: A classical molecular dynamics study[J]. AIP Conf Proc, 2012, 1 496 (1): 229-232.

[14] 唐杜. 硅基器件的单粒子翻转及单粒子位移损伤的数值模拟研究[D]. 西安：西安交通大学, 2016. (TANG Du. Numerical simulations on single event upset and single particle displacement damage in silicon devices[D]. Xi’an: Xi’an Jiaotong University, 2016.)

[15] FICHTHORN K A, WEINBERG W H. Theoretical foundations of dynamical Monte Carlo simulations[J]. J Chem Phys, 1991, 95 (2): 734-750.

[16] MARTIN-BRAGADO I, ZOGRAPHOS N. Indirect boron diffusion in amorphous silicon modeled by kinetic Monte Carlo[J]. Solid-State Electronics, 2011, 55 (1): 25-28.

[17] MARTIN-BRAGADO I, TIAN S, JOHNSON M, et al. Modeling charged defects, dopant diffusion and activation mechanisms for TCAD simulations using kinetic Monte Carlo[J]. Nucl Instrum Methods Phys Res B, 2006, 253 (1/2): 63-67.

[18] FISICARO G, PELAZ L, ABOY M, et al. Dopant dynamics and defects evolution in implanted silicon under laser irradiations: A coupled continuum and kinetic Monte Carlo approach[C]// International Conference on Simulation of Semiconductor Processes and Devices (SISPAD). Glasgow, Scotland, 2013: 33-36.

[19] TANG D, MARTIN-BRAGADO I, HE C H, et al. Time dependent modeling of single particle displacement damage in silicon devices[J]. Microelectronics Reliability, 2016, 60: 25-32.

[20] 唐杜, 贺朝会, 臧航, 等. 硅单粒子位移损伤多尺度模拟研究[J]. 物理学报, 2016, 65(8): 084209. (TANG Du, HE Chao-hui, ZANG Hang, et al. Multi-scale simulations of single particle displacement damage in silicon[J]. Acta Phys Sin, 2016, 65(8): 084209.)