0 引言

弥散在石墨集体中的TRISO(Tristructural-Isotropic)包覆颗粒燃料是高温气冷反应堆的基本元件, 从内到外它的基本组成单元依次是:燃料内核(通常是UO$_{2})$、多孔石墨缓冲层、内致密热解炭层、SiC层、外热解炭层, 如图 1所示.在TRISO颗粒中SiC层为其提供了结构支持并且起到阻碍主要核裂变产物扩散的作用^[1].虽然SiC层能够有效抑制大部分金属裂变产物的扩散, 但是也有实验发现, 在高温条件下仍然有Ag等裂变产物透过SiC层扩散到反应堆冷却剂中^[1-2].这会导致Ag在反应堆中的积淀, 这给核电系统的维护工作造成了很大的困难, 同时还存在一定的安全隐患, 因为放射性产物$^{110}$Ag的半衰期比较长, 为253d.

为了控制裂变产物Ag在SiC层中的扩散, 明确其扩散机制是至关重要的.针对Ag在SiC层中的扩散现象, 很多扩散机制被提出, 其中包括SiC的降解^[3-4]、通过裂缝和纳米孔的扩散^[5-7]、沿着晶粒边界和错位层的扩散^[8-9]等.但是时至今日, 究竟是哪一种扩散机制起作用仍然是不明确的, 有待于实验和理论工作的进一步验证.

关于Ag在SiC中的扩散, 在实验^[10-14]和基于密度泛函理论(Density Functional Theory, DFT)的第一性原理^[15-16]方面都有相应的研究, 其中对Ag扩散系数的计算都是根据计算出来的主要几种点缺陷的形成能和迁移势垒按照公式$D=D_0\exp(-Q/k_bT)$估算获得; 然而没有分子动力学模拟方面的研究, 这是因为缺乏相应的原子间相互作用势.众所周知, 由于计算条件的限制, 第一性原理方法只能在小尺度范围内进行计算研究; 而分子动力学方法可以在大尺度下进行长时间的动态模拟, 并且可以通过分子动力学模拟对扩散机制的微观过程进行直接观察.分子动力学方法可以用均方位移法或速度关联函数法直接计算扩散系数.对于缺陷迁移这类微观上低概率事件, 可以借助于加速分子动力学方法^[17]来进行模拟.

进行分子动力学模拟研究需要相应的势函数作为支持, 在AgSiC三原子体系中存在Si-Si、C-C、Si-C、Ag-Ag、Ag-Si和Ag-C等6种相互作用关系.但是, 现有的势中没有Ag-Si和Ag-C势.况且, 已有的Si-Si势、C-C势、Si-C势^[18-22]和Ag-Ag势^[23]并不是针对我们研究的体系拟合的, 需用进一步修正甚至重新拟合.因此, 需要建立一个AgSiC三原子体系的相互作用势, 为分子动力学方法研究Ag在SiC中扩散提供支持.

1 理论方法 1.1 计算原理及公式

本文中势函数的验证过程中所涉及的体系内聚能、缺陷形成能和体弹性模量等的计算分别采用如下公式.

内聚能的计算公式^[24]为

$E_\mathrm {coh}=E_\mathrm {bulk}-\sum\limits_in_iE_\mathrm {atomi},$

(1)

其中, $E_\mathrm {bulk}$和$E_\mathrm{atomi}$分别表示单个化学势晶体能量和单个原子能量(其中包括自旋能和零点振动能两部分).

缺陷形成能的计算公式^{[15, 17]}为

$E_f=E_\mathrm {def}-E_\mathrm {undef}+\sum\limits_I\Delta n_I\mu_I, \quad\mu_I=(E_I+\gamma_I),$

(2)

其中, $E_\mathrm {def}$表示包含缺陷的晶胞能量, $E_\mathrm{undef}$表示完整晶胞能量, $I$表示元素种类, $\Delta n_I$表示完整晶胞中$I$类原子数减去缺陷晶胞中$I$类原子数, $\mu_I$表示对应化合物构型(SiC晶体)下$I$类原子的化学势, $E_I$是相应单质构型(在这里就是Ag、Si和C晶体的)下单个$I$原子的能量, $\gamma_I$表示元素$I$在相应单质构型下的化学势, 可根据文献[25]确定它的取值范围, 即, $0\leq\gamma_I\leq E_\mathrm{SiC}-E_\mathrm{Si}-E_\mathrm C$.这里分3种情况分别取^[15]:$\gamma_\mathrm{Si}=\gamma_\mathrm C=(E_\mathrm{SiC}-E_\mathrm{Si}-E_\mathrm C)/2$; $\gamma_\mathrm{Si}=E_\mathrm{SiC}-E_\mathrm{Si}-E_\mathrm C, \gamma_\mathrm C=0; \gamma_\mathrm{Si}=0, \gamma_\mathrm C=E_\mathrm{SiC}-E_\mathrm{Si}-{E_\mathrm C}$.本文计算SiC晶体缺陷形成能时只考虑γ_C=(E_SiC−E_Si−E_C ^[15])/2的情况.

晶格常数$a$、内聚能E、体弹性模量$B_{0}$等值是通过Murnaghan方程^[26]对体积-能量关系的拟合结果获得, Murnaghan方程的一般表达式为

$E(V)=\frac{B_0V}{B_0'(B_0'-1) }\Bigg(B_0'\Bigg(1-\frac{V_0}{V}\Bigg)+\Bigg(\frac{V_0}{V}\Bigg)^{B_0'}-1\Bigg)+E(V_0) .$

(3)

本文采用Tersoff势^[21]来描述原子间相互作用, 其一般形式为

$E=\Sigma_iE_i=\frac{1}{2}\Sigma_{i\neq j}V_{ij}, \quad V_{ij}=f_C(r_{ij})(f_R(r_{ij})+b_{ij}f_A(r_{ij})), $

(4)

其中, 排斥(Repulsive)势$f_{R}$和吸引(Attractive)势$f_{A}$如式(5) 所示, 截断(Cutoff)函数如式(6) 所示, 即

$ f_R(r_{ij})=A_{ij}\exp(-λ_{ij}r_{ij}), \quad f_A(r_{ij})=B_{ij}\exp(-\mu_{ij}r_{ij}), $

(5)

$ {f_C}({r_{ij}}) = \left\{ {\begin{array}{*{20}{l}} {1,\quad {r_{ij}} ＜ {R_{ij}},}\\ {\frac{1}{2} + \frac{1}{2}\cos [\pi ({r_{ij}} - {R_{ij}})/({S_{ij}} ＜ {R_{ij}})],{R_{ij}} ＜ {r_{ij}} ＜ {S_{ij}},}\\ {0,\quad {r_{ij}} > {S_{ij}},} \end{array}} \right. $

(6)

以及表示原子周围环境的三体作用项

$\begin{align} & b_{ij}=\chi_{ij}(1+\beta_I^{n_i}\zeta_{ij}^{n_i})^{-1/2n_i},\quad \zeta_{ij}=\sum\limits_{k\neq i,j}f_C(r_{ik})g(\theta_{ijk}),\notag\\ & g(\theta_{ijk})=1+c_i^2/d_i^2-c_i^2/(d_i^2+(h_i-\cos\theta_{ijk})^2). \end{align}$

(7)

式(4) 至式(7) 中$i, j, k$表示原子序号, $r_{ik}$表示原子i和k之间的距离, $\theta_{ijk}$表示$ij$和$ik$原子间化学键之间的角度.为了简单起见, 令$\chi=1$, 则以上各式共有11个待定参数.其中$A$, $B$, $λ$和$\mu$这4个参数为两体相互作用部分, $\beta, n, c, d$和$h$为三体相互作用部分, $R, S$的取值跟截断半径有关.

在拟合过程中POTFIT采用最小二乘法、共轭梯度法、模拟退火法等, 来优化参考构型的力和能量等物理量的第一性原理计算值和势函数计算值之间的差值函数, 使之取最小值, 从而获得势函数各参量的值.差值函数的一般形式为

$Z=Z_E+Z_F+Z_S, $

(8)

其中,

$\begin{align*} & Z_E=\sum\limits_{i=1}^{N_C}W_i\frac{(E_i^\mathrm {tersoff}-E_i^\mathrm{DFT})^2}{(E_i^\mathrm{DFT})^2}, \\[2mm] & Z_F=\sum\limits_{i=1}^{N_A}\sum\limits_{j=1}^{3}W_i\frac{(F_{i, x_j}^\mathrm{tersoff}-F_{i, x_j}^\mathrm {DFT})^2}{(F_{i, x_j}^\mathrm{DFT})^2}, \\[2mm] & Z_S=\sum\limits_{i=1}^{N_C}\sum\limits_{j=1}^{6}W_i\frac{(\sigma_{i, j}^\mathrm{tersoff}-\sigma_{i, j}^\mathrm{DFT})^2}{(\sigma_{i, j}^\mathrm{DFT})^2}, \end{align*}$

分别代表来自内聚能、力和压强的贡献.在以上各式中$W_i$表示权重, $N_C$表示所有构型的数目, $N_A$表示所有原子的数目, $F_{i, x_j}$表示i原子在$x_j$方向受到的力, $\sigma_{i, j}$表示i原子应力张量的第j个数值.本文拟合时只考虑了式(8) 中的能量和力两项对差值函数的影响.

1.2 计算参数设置

本文采用Tersoff势作为Ag、Si和C原子间作用的力场函数, 其合理性已经在参考文献[18-22, 27]中得到了很好验证.拟合势函数所需的第一性分子动力学、静态优化和电子步优化等计算及验证势函数所需的弹性常数、缺陷形成能、晶格常数和内聚能等第一性原理的计算采用基于密度泛函理论(DFT)的VASP^[28]软件包.势函数的拟合采用基于“力匹配"方法的POTFIT^[29]软件包.势函数验证过程中用拟合的势计算弹性常数、晶格常数、内聚能、体弹性模量和SiC晶体缺陷形成能等值时采用基于“牛顿力学"的Lammps软件包.

本文用VASP软件包的投影缀加波方法进行密度泛函计算, 采用PBE^[30]形式的GGA^[31]交换关联泛函求解Kohn-Sham方程.平面波截断能取570 eV. Monkhorst-Pack K-point^[32]网格:8原子体系取13$× $13$× $13;大于60原子体系取7$×$7$× $7.在SiC和AgSiC缺陷的第一性分子动力学计算中, 步长取1.55 fs, 温度参考文献[10-14]中的实验温度取1073~1 773 K.

Lammps计算过程中涉及的静态优化过程都采用最陡下降算法, 步长取1.55fs. POTFIT优化过程中, 体系能量权重取100, 退火温度取1 500 K.

2 结果与讨论 2.1 构型选取

在势函数的拟合过程中要选取包含所研究的物理过程或物理条件的晶体结构作为第一性原理计算的构型.为了描述Ag在SiC晶体中的扩散, AgSiC三原子缺陷体系构型的选取David Shrader在文献[15]中提到的16种缺陷作为基本构型, 它们分别是Ag$_\mathrm{C}$、Ag$T_\mathrm{C}$、Ag$_\mathrm{C}V_\mathrm{C}$、Ag$_\mathrm{C}V_\mathrm{Si}$、Ag$_\mathrm{C}2V_\mathrm{Si}$、Ag$_\mathrm{C}V_\mathrm{Si}V_\mathrm{C}$、Ag$_\mathrm{C}V_\mathrm{Si}$Si$_\mathrm{C}$、Ag$_\mathrm{C}V_\mathrm{Si}$C$_\mathrm{Si}$、Ag$_\mathrm{Si}$、Ag$T_\mathrm{Si}$、Ag$_\mathrm{Si}V_\mathrm{Si}$、Ag$_\mathrm{Si}V_\mathrm{C}$、Ag$_\mathrm{Si}2V_\mathrm{C}$、Ag$_\mathrm{Si}V_\mathrm{C}V_\mathrm{Si}$、Ag$_\mathrm{Si}V_\mathrm{C}$Si$_\mathrm{C}$和Ag$_\mathrm{Si}V_\mathrm{C}$C$_\mathrm{Si}$.在这里Ag$_\mathrm{C}$和Ag$_\mathrm{Si}$分别表示Ag原子替换了C和Si原子, $V_\mathrm{Si}$和$V_\mathrm{C}$分别表示在Ag原子最近的Si空位和C空位, 它们前面的数字表示空位的数目; Si$_\mathrm C$和C$_\mathrm {Si}$则分别表示1个Si原子替换了离Ag原子最近的C原子和C原子替换了离Ag原子最近的Si原子, Ag$T_{\mathrm C}$和Ag$T_\mathrm{Si}$则分别表示Ag填隙在4个C的中心和4个Si的中心.为了便于理解, 我们画出了其中5种缺陷构型, 如图 2所示, 其中, 图 2(a)Ag$_\mathrm {Si}V_\mathrm{C}$表示Ag原子取代了Si原子的位置, 并且在Ag原子的4个最近邻位置之一出现了1个C空位; 图 2(b)Ag_SiV_CV_Si则表示取代Si的Ag原子的四个最近邻位置上出现1个C空位, 同时在Ag次近邻位置上出现了一个Si空位; 图 2(c)Ag$_\mathrm{Si}2V_\mathrm{C}$表示取代Si的Ag原子的4个最近邻位置上出现2个C空位; 图 2(d)Ag$_\mathrm{Si}V_\mathrm{C}$C$_\mathrm{Si}$表示取代Si的Ag最近邻位置上出现1个C空位, 同时次近邻位置的Si被C取代; 图 2(e)Ag$_\mathrm{Si}V_\mathrm{C}$Si$_\mathrm{C}$则表示取代Si的Ag最近邻位置上出现1个C空位, 同时最近邻位置的另一个C被Si原子取代.其他11种构型表示方法以此类推.

本文研究的AgSiC三原子体系的主体是SiC材料, 所以在势的拟合过程中要兼顾SiC晶体的一些性质.因此我们首先对Si-Si, C-C和Si-C原子间作用势进行了拟合.在本文中Si-Si、C-C和Si-C原子间Tersoff势的截断半径^[21]分别取3.0 Å, 2.1 Å和2.5 Å; 而在完全优化后的SiC晶格中Si-Si, C-C和Si-C原子间距分别约为3.1 Å、3.1 Å和1.9 Å.Si-Si和C-C间距都为3.1 Å大于对应的势函数截断半径3.0Å和2.1 Å, 因此在计算完美SiC晶体性质时, Si-Si和C-C势对其结果几乎没有影响; 而Si-C间距1.9Å在Si-C势的截断半径2.5 Å之内, 因此Si-C作用在SiC晶体性质计算时有很大影响.我们用文献[21]中给出的Tersoff势进行了验证, 结果表明, 只有在SiC晶体中存在间隙等缺陷时, Si-Si和C-C作用势才会对SiC晶体性质的计算结果产生影响, 否则没有影响.这是因为间隙的存在使得Si-Si和C-C原子间距分别在3.0Å和2.1 Å以下产生了分布.因此, 我们认为Si-Si和C-C作用势对SiC晶体性质直接产生的影响很小, 可以最先单独进行拟合, 这样使得拟合过程快速和简单, 并且得到的SiC体系的势还可以兼顾C和Si晶体各自的性质.

基于以上分析, 具体拟合过程分3步进行:首先, 分别拟合Si-Si和C-C势; 然后, 固定Si-Si和C-C势的各参数不变进一步拟合Si-C势; 最后, 固定前两步中的拟合结果不变, 进行Ag-Si和Ag-C作用势的拟合.相应构型的选取也分为如下3类.

1. Si-Si和C-C作用势构型

(1) 首先, 完全优化64原子Si(C)晶格; 固定晶格常数分别对65(间隙)、63(空位)原子的Si(C)构型进行等体积优化, 并计算相应构型Si(C)原子间距的分布及所对应权重(这里的权重与相同间距出现的次数成正比). Si和C单质体系的间隙、空位结构如图 3所示.

(2) 根据(1) 中原子间距分布, 放大或缩小8原子完全优化的晶格, 使其原子间距刚好满足(1) 中得到的原子间距分布.一般会产生10至20个构型.

2. Si-C作用势构型

(1) 状态方程.为了确定SiC晶体的晶格常数、形成能、体弹性模量和弹性常数, 我们对SiC完整晶格进行了完全优化(平衡体积为$V_{0})$, 并选取了0.90$V_{0}$至1.15$V_{0}$的15个构型.

(2) 缺陷.为了初步考虑AgSiC体系中的缺陷, 我们对Ag$_\mathrm{C}$等16种缺陷固定体积进行了优化(采用64原子SiC晶体完全优化的晶格常数), 然后去掉缺陷构型中Ag原子并保持原有Si原子、C原子位置不变, 进行DFT静态计算(只做电子步计算).缺陷结构参照图 2.此外, 我们还考虑了SiC中的本征空位$V_{\mathrm C}$和$V_\mathrm{Si}$两种构型, 其结构如图 4所示.

(3) 分子动力学过程.为了获得SiC晶体在相应温度下可能出现的构型及其力和能量等信息, 我们取64原子的SiC超胞, 采用NVT系综在1 073~1 773K下进行第一性原理分子动力学计算, 每个分子动力学过程进行500步, 步长1.55 fs, 从轨迹中选取20个不同温度下的构型.

3. Si-Ag和C-Ag作用势构型

(1) 缺陷静态优化.对Ag$_\mathrm {C}$等16种缺陷进行固定体积优化, 其晶格常数及体积取完全优化的SiC晶体的值.

(2) 缺陷分子动力学过程.为了获得各缺陷体系在相应温度下可能出现的构型及其力和能量等信息, 我们用第一性原理分子动力学方法对Ag$_\mathrm{C}$等16种缺陷在NVT系综下运行500步, 步长1.55 fs, 温度1 073$\sim $1 773 K.每种缺陷在不同温度下选取10个构型.

2.2 拟合过程

拟合过程中, 首先, 根据第2.1节中所述的3类构型选取构型, 并用第一性原理计算各构型的力和能量等信息.然后, 用POTFIT最优化Tersoff势函数中的11个待定参数, 调整各参数的变化范围直到结果收敛.最后, 用所得的Tersoff势计算各体系晶格常数$a$、内聚能E、弹性常数$C_{ij}$、体弹性模量$B_{0}$、缺陷形成能和16种缺陷的内聚能, 并跟第一性原理计算结果进行对比, 根据对比结果调整各个构型的权重, 反复拟合, 直到对比结果相吻合.

在Si-Si和C-C势的拟合过程中, 我们只对B、$λ$和$\mu $这3个参数进行了改动, 其他参数借鉴文献[21]中的数值保持不变.之所以这么做是由于本文拟合基于第一性原理计算, 而文献[21]中的拟合基于实验数据, 两者能量不匹配, 而$A$、$B$、$λ $和$\mu $等参数主要确定体系能量的高低, 这里$B$、$λ$和$\mu $参数足以达到能量匹配的效果; $R$和S参数确定相互作用势的截断范围, $\beta $、n、c、$d$和$h$跟原子周围环境因素有关, 都可以不用改变.

在拟合Si-C势的拟合过程中, 除了表示截断的$R$和S这2个参数与文献[21]保持一致外, 其他参数都进行了改变.这是由于除了能量不匹配外, 本文考虑的缺陷与文献[21]中考虑的缺陷不完全同, 从而部分原子的周围环境因素不完全相同, 跟周围环境因素有关的参数需用重新拟合.

在Ag-Si和Ag-C作用势的拟合过程中, $\beta$和n参数的选取参考了文献[33]中的结果直接取1, $R$和S的取值结合SiC晶体中原子最近邻和次近邻间距进行调整.

2.3 势的验证及结果

由第一性原理计算的Si、C和SiC晶体的晶格常数、内聚能、体弹性模量、弹性常数、缺陷形成能等物理量见表 1DFT所示.这些数据作为势函数验证的主要依据, 其准确性决定验证环节的可靠程度.因此我们拿这些数据跟其他第一性原理计算结果(表 1Refa和Refb)进行了对比, 由表 1DFT、Refa和Refb的数据可以看出, 本文的第一性原理计算值跟其他第一性理论计算值非常接近, 其误差跟计算条件及精度有关, 可以接受.

图 5中(a)、(b)和(c)分别给出了Si、C和SiC晶体的内聚能-体积关系及用Murnaghan方程拟合的结果.图中, 纵坐标为平均每个原子的内聚能, 横坐标为平均每个原子的体积, 黑色“$\triangle$"、蓝色“$\Box$"和红色“$\bigcirc$”分别表示文献[21]中所给出的势(Tersoff)、第一性原理(DFT)和本文拟合的势(This work)的计算值, 相应的同种颜色曲线代表Murnaghan方程拟合结果.由图 5可以看出, 由第一性原理数据的Murnaghan曲线(蓝色)和本文给出的势的数据的Murnaghan曲线(黑色)基本重合, 两者符合得很好.由Murnaghan方程拟合得到的晶格常数、内聚能和体弹性模量的值见表 1所示.

表 1给出了用本文拟合得到的势计算获得的Si、C和SiC晶体中晶格常数$a$、内聚能E、体弹性模量$B_{0}$、弹性常数$C_{ij}$、空位形成能Vac、间隙形成能Int、Si空位形成能$V_\mathrm{Si}$和C空位形成能$V_{C}$等物理量数值, 并跟第一性原理计算值、实验值及别人第一性原理计算值等对比.其中, exp、Tersoff和Ref分别表示实验值、文献[21]中势函数的计算值和参考文献的第一性原理计算值, 相应的数据来源已在表格中标出; This work和DFT分别表示本文拟合所得势的计算值和第一性原理计算值; 空位Vac、间隙Int、Si空位$V_\mathrm{Si}$和C空位$V_\mathrm {C}$的结构如图 3和图 4所示.由表 1(a)和表 1(b)可以看出, 本文拟合的势对Si和C晶体的晶格常数、内聚能、弹性常数和体弹性模量的计算值跟第一性原理计算值非常吻合; 而对Si和C晶体中空位和间隙形成能的计算值跟文献[21]给出势的计算值比较接近, 但与其相比更接近第一性原理计算值, 其结果有所改善; 这是由于我们在Si-Si和C-C势的拟合过程中只拟合了$B$、$λ$和$\mu $这3个参数, 同时在构型选取时考虑了缺陷体系的原子间距分布的缘故.由表 1(c)SiC可以看出, 本文拟合的势对其他物理量的计算值跟第一性原理计算值吻合得很好.

表 2给出了计算AgSiC缺陷形成能所需的化学势的值, 并与文献[15]中的数据进行了对比, 其计算方法及物理意义已在第1.1节中介绍.由第1.1节中的介绍可知, 本文和文献[15]中给出的各化学势的差异归根结底是$E_\mathrm{Si}$、$E_\mathrm {C}$和$E_\mathrm {SiC}$值的不同.文献[15]中取块体能量-5.44 eV、-9.20 eV和-15.08eV作为它们的值, 本文取内聚能-4.545 eV、-7.773 eV和-12.862 eV作为它们的值.考虑Si和C的原子能-0.880eV和-1.316 eV后, 根据式(1) 可由内聚能计算块体能, 其值分别为-5.425 eV、-9.089 eV和-15.059 eV.这跟文献[15]中给出的值仍然有微小差别, 这是由于计算精度及方法不同所致.

表 3给出了本文拟合的势对16种AgSiC三原子体系缺陷形成能的计算值跟第一性原理计算值及文献[15]中给出的值的对比结果及相对误差.可以看出, 本文拟合的势计算16种AgSiC体系缺陷形成能的值与第一性原理计算的值符合得很好.除了Ag$_\mathrm CV_\mathrm {Si}V_\mathrm C$、Ag$_\mathrm CV_\mathrm {Si}$Si$_\mathrm C$、Ag$_\mathrm CV_\mathrm C$和Ag$_\mathrm {Si}V_\mathrm C$Si$_\mathrm C$4这种缺陷形成能误差分别为2.12 eV、2.31 eV、4.83 eV和2.31 eV外, 其他缺陷形成能的误差都不超过0.84 eV.由于我们考虑的缺陷种类较多, 缺陷结构差异较大, 拟合得到的势很难顾及到每一个缺陷的性质, 这是经验势本身的局限性所致.个别构型的缺陷形成能显得大, 但这也是可以接受的.如参考文献[20-21]中的势对int(S)和$V_\mathrm {Si}+V_\mathrm{C}$缺陷形成的计算值跟第一性原理计算值的最大误差分别为6.7eV和5.3 eV, 相对误差分别为40%和42%.

最终得到的Si、C和Ag原子间相互作用势参数如表 4所示, 本文所研究的缺陷体系中只有一个Ag原子, 故暂不考虑Ag原子间相互作用, Ag-Ag势参数手动设置使之为零.我们用表 4中的势分别计算了Si、C和SiC的多个物理量并跟实验值(exp)进行了对比, 如表 1所示, 两者有一定的差别, 这是由于我们的拟合是基于第一性原理计算进行的.这个误差的根源是第一性原理计算结果和实验结果之间的差别.总体上来说我们的拟合是成功的.

3 结论

本文采用“力匹配”方法, 基于第一性原理计算拟合了存在辐照缺陷的SiC晶体结构中Ag、Si和C原子间相互作用势.这种方法避免了传统的基于实验数据拟合势的方法中数据源不足的问题, 为复杂体系材料原子间作用势的拟合提供了可能.拟合结果表明, 通过该方法获得的原子间相互作用势能够很好地描述Si、C和SiC晶体材料的晶格常数、内聚能、弹性模量和缺陷形成能等.采用Tersoff形式的势来描述Ag掺杂的SiC晶体结构中Ag、Si和C原子间相互作用行之有效, 不但可以兼顾Si、C和SiC晶体各自的性质, 而且还能准确计算AgSiC体系的缺陷形成能.当然, 由于经验势本身的局限性, 我们对缺陷形成能的拟合中, 个别构型的误差还是较大, 这方面还有一定的提升空间.为拟合更为精确的势, 需要更加全面地去考虑各个因素的影响, 我们的现有结果为进一步的研究奠定了初步基础.