高精度拟合（LAMMPS计算MSD拟合原子扩散能垒示例）

跟据过度态理论，体系中某一group原子的扩散系数与均方位移具有以下公式（1）和（2）的关系，公式（1）中r(t)为原子位移，d为维度，t为测量原子扩散的时间。公式（3）为扩散系数与温度T的关系，其中E_a为扩散能垒，k_B为波尔滋蔓常数。因此只需要计算一定温度范围内扩散系数与温度的关系，就可以通过公式（3）拟合扩散能垒。

本文笔者将利用LAMMPS以H在Si@SiO2异质结处扩散为例，示范通过过渡态求300K，H扩散系数的过程。同时笔者会给出in文件，读者可经过简单的修改，求出一定温度范围的扩散系数与温度的数值关系从而拟合扩散能垒。

in文件脚本如下：

### Basic setup ###variable t equal 300units metaldimension 3boundary p p patom_style atomicnewton on### Structure ###read_data pos1.datneighbor 1 bincreate_atoms 3 single 5 5 22### Potentials ###pair_style hybrid tersoff lj/cut 10.0pair_coeff * * tersoff SiO.tersoff Si O pair_coeff 1 3 lj/cut 0.0057673 3.1988pair_coeff 2 3 lj/cut 0.0022281 2.8446pair_coeff 3 3 lj/cut 0.001908 2.5711mass 1 28.085mass 2 16mass 3 1.008### Group ###group H type 3### Minimization ###thermo 100thermo_style custom step time vol atoms temp ke pe lx pressfix 1 all box/relax iso 0.0 vmax 0.01min_style cgminimize 1e-25 1e-25 100000 100000### NPT1 ###reset_timestep 0thermo 1000thermo_style custom step dt time atoms temp ke pe vol press lx ly lz velocity all create $t 92137453 dist gaussianfix NPT all npt temp $t $t 1 iso 0.0 0.0 1.0 drag 0.2timestep 0.001run 10000reset_timestep 0compute Msd H msd fix s all vector 1 c_Msd[4]variable MSD equal sum(f_s)/1000 variable stp equal step*dt### NPT2 ###fix msd1 all print 1000 "${stp} ${MSD}" file msd_$t.txt screen nothermo 1000thermo_style custom step temp ke pe vol press lx v_MSDtimestep 0.001run 500000unfix NPT

运行完毕后通过对计算所得MSD与时间得关系进msd_300.txt文件拖入Oringlab行拟合，可得300 K下H在Si@SiO2异质结处扩散系数为8.295 Ų/ps=8.295*10⁸m²/s。

--------------------------------------------------------------

PS：由于时间仓促，可能有错误之处，欢迎各位同行交流学习，同时若有疑问、计算技术讨论、合作可以联系笔者邮箱：1924311399@qq.com谢谢！

公众号推荐：计算运维鸟