锂电池材料Li2O力场开发及性质研究
—力学性质、热性质、输运性质
1. 研究背景
当今，锂电池作为新能源材料在很多领域中都具有极为广泛的应用，比如微电子设备、电动汽车、航空航天、军事设备等。为了加快新电池材料的开发设计周期，必须要理论与实验相结合，运用可靠的计算模拟工具实现多尺度理论研究支撑实验工作。Li2O是锂离子电池中的一种基本氧化物材料，因此，对Li2O中Li离子的扩散行为、热性能和输运性质等的研究极为重要。但是，在关于Li2O的分子动力学研究中，仍旧缺少可靠的基础力场参数，而在此案例中，先确定了适合于Li2O的力场参数，然后利用此力场对Li2O进行了理论计算研究。
在本案例中，开发MedeA软件的原厂美国Materials Design 公司与日本Toyota Central Research and Development Laboratories从理论和实验两个角度联合研究了Li2O这一新电池材料的特性。其中，作者首先采用MedeA–VASP模块中的从头分子动力学确定了Li2O的力场参数，开发了适用于用Li2O电池材料的力场；然后，采用新开发的力场AICl3，结合MedeA–LAMMPS、MedeA–Mechanical Thermal（MT）、MedeA–LAMMPS-Diffusion模块对此体系的机械性能和扩散性质做了系统地研究。研究表明，AICl3力场能够通过设置较少参数项就达到预测热性能的目的，且具有足够的准确度。具体结果如下：
2. AICl3力场参数计算
作者采用MedeA-VASP模块中的GGA-PBEsol泛函优化Li2O晶胞，然后通过Supercell Builder超晶胞建模工具创建(Li2O)16 超晶胞。接着，采用MedeA-VASP模块中的从头分子动力学对超晶胞模型进行NVE系综计算。由NVE计算得到的总能需要进行偏移校对，而此总能偏移量又与重要力场参数Li离子的有效电荷具有一定的对应关系。
图1 计算的熔点随总能偏移量变化关系 图2 Li离子有效电荷随总能偏移量变化关系
图一为计算得到的熔点与总能偏移量之间的关系，得到线性关系后作者根据实验中确定的Li2O熔点确定最终总能偏移量（-2.47 eV/atom）。图二为总能偏移量与Li离子有效电荷之间的线性关系，作者由总能偏移量确定Li离子的有效电荷（0.79｜e︱）。Li离子有效电荷的改变将对整个体系的排斥作用项参数产生影响。
表一为经过熔点校正的力场（AICl3）和FIT-GGA（由Oda等人的实验数据确定）力场参数对比。采用AICl3力场NPT系综下的模拟结果如下：较低温度显示为晶体状态，较高温度显示为熔化状态。
表 1 Li2O力场参数（通过熔点校正）
3. AIC13力场测试分析
图3显示了AICl3力场中Li、O原子间的相互作用。Li-O排斥作用的范围在2 Å 以下，这与Li2O固体中Li-O平衡的距离 1.99 Å十分吻合。
图4中实线为VASP从头分子动力学计算结果，虚线为基于分子力学采用AICl3力场的计算结果，可以发现两者十分一致。这说明AICl3力场中的有效电荷与VASP计算得到的Bader Charge也具有非常高的吻合度。图5表明由AICl3力场计算的结果比由FIT-GGA力场计算的结果与VASP的计算结果更接近。
图 3 由AICl3 力场计算的Li离子与O离子相互作用势
图 4 (Li2O)16超晶胞于1000 K下分别采用VASP和AICl3计算得到的能量超曲面
图 5 从(Li2O)16超晶胞中移走一个Li离子，将体系中其他原子均固定，计算其于3000K下Li离子空穴附近的能量超曲面。分别采用VASP、FIT-GGA和AICl3力场计算
由图3、图4和图5的结果表明，作者新开发的AICl3力场对于Li2O电池材料体系更合适。
4. 由AICl3力场计算的Li2O理化性质
4.1 晶格常数和力学性质
作者通过四种不同的方法计算了Li2O的晶格常数，包括量子力学（MedeA-VASP模块）和分子力学（MedeA-LAMMPS模块采用FIT-GGA、AIM和AICl3三种力场）；通过MedeA-MT模块计算了体系的力学性质。具体结果如表2。结果表明，计算得到的晶格常数都比实验值小；密度比实验值大；对于弹性常数和弹性模量，AICl3与FIT-GGA的计算结果相似；但对于熔点，只有AICl3力场的计算结果在误差范围内。
表 2 Li2O 材料实验及计算性质。结构参数由MedeA-VASP计算。弹性系数、熔点由MedeA-LAMMPS中的FIT-GGA力场计算
4.2 热膨胀
接着，作者采用MedeA-LAMMPS模块基于AICl3力场计算了Li2O材料的热膨胀行为。
图6呈现了密度描述为温度的函数，很好的说明了此材料的热膨胀行为。在T = 0 K时， 由AIC13力场计算的密度介于第一性PBEsol泛函和LDA泛函结果之间；在Tc = 1200 K时，密度明显减小，这是由于次晶格位的Li离子发生预熔（pre-melting）现象造成的，此时，Li离子液会呈现一定的无序性；在Tmelt = 1711K时，材料密度减小了约 12 %，并随着温度的升高以线性关系下降。
图 6 Li2O的密度同温度的对应关系图。实心点表示了采用三种不同的DFT泛函。其中，PBEsol泛函最接近0 K下的实验值
图 7 由AICl3力场（虚线）计算得到的热膨胀曲线与实验值（符号）的对比图。
图7为通过AICl3力场计算的热膨胀∆L/L0曲线同实验值的对比。结果表明，AICl3力场同实验值具有很好的拟合度。在T = 700 K时，计算得到的热膨胀效果比实验平均值略高；在T > 1100 K时，实验值与计算值的斜率都有明显增加。
图8显示了热膨胀系数关于温度的函数曲线，同样为基于AICl3力场计算值与实验值的对比结果。从低温到900 K的升温过程中，热膨胀系数随着温度的升高而升高，与实验值充分吻合。在1711 K熔点附近，热膨胀系数的斜率与 900 K到1200 K的斜率相同，而比900 K以下的斜率要大，说明在较高温度的液相中，热膨胀的效果更强。
图 8 计算与实验的热膨胀系数关于温度的函数。（1）、（2）分别为实验值。T = 1200 K的峰体现了次晶格位置的Li离子在超离子导体中的预熔现象。
4.3 扩散
然后，作者又通过MedeA-LAMMPS-Diffusion研究了三种不同Li2O体系的扩散系数关于温度的函数变化曲线：无空穴、5×5×5的 (Li2O)500超晶胞；含有0.4 %Li空穴和0.2 %O空穴的超晶胞模型(Li2O)499；以及含有2 %Li空穴和1 %O空穴的超晶胞模型 (Li2O)495。
图 9 扩散系数关于温度的函数曲线（计算和实验结果对比）
图9表明，当温度升高时，Li离子的扩散速率加快。同Oda开发的FIT-GGA力场计算结果相比，AICl3计算得到的扩散系数更接近实验结果。此外，随着Li离子空穴浓度的增大，其扩散系数也更大，这与Oda的计算结果和实验结果完全一致。当温度达到1100 K以上时，体系中空穴对其扩散行为的影响减小，这与次晶格位Li离子的预熔现象有关。
为了更加清楚地了解的次晶格位Li离子的预熔问题，作者又进一步分别分析了Li离子和O离子关于温度的扩散行为，计算了不同温度下的均方位移。如图10所示。
图 10 Li2O中Li离子和O离子的均方根位移关于时间的函数变化 。可以看到在900- 1000 K之间Li离子开始扩散，并随着Li2O的液态化，其Li离子的扩散程度增加。
通过Li和O离子的MSD均方根位移的变化可知，Li离子在900- 1000 K时开始扩散，并在1200 K时达到热扩散系数的峰值，同图8的结果完全吻合，这充分说明了液态Li的形成。不同的是，O离子始终保持在晶格原来的位置，即使到了熔点1711 K也不发生变化。
5. 结论
本案例中，Materials Design公司通过第一性计算数据的校正，开发了适用于研究Li2O材料的新力场AIC13。经过与Toyata开发实验室获得的实验数据对比，AIC13在计算Li2O热性质上，比其他力场与实验值具有更高的匹配度。同时，作者采用MedeA中多个模块基于AICl3力场对Li2O电池材料的力学、热膨胀、扩散、融化等性质进行了系统研究，均得到了与实验值十分吻合的结果。因此，通过MedeA计算及性质预测模块的结合，能够帮助预测锂电池材料的多种热性质，从而缩短其开发设计周期。
参考文献：
R. Asahi, C. M. Freeman, P. Saxe, E. Wimmer, Thermal expansion, diffusion and melting of Li2O using a compact forcefield derived from ab initio molecular dynamics, Modeling and Simulation in Materials Science and Engineering, 2014, 22: 075009
使用MedeA模块：
Welcome to MedeA Bundle
MedeA - VASP
MedeA - LAMMPS
MedeA - Forcefields
MedeA - MT
MedeA - LAMMPS-Diffusion
