电池材料 LiFePO4 中锂离子的扩散

<chem>LiFePO4</chem> 是典型的可重复充电的锂离子电池电极材料，成本低，环境友好^[1,3]。在本教程中，将使用 ATK-DFT 估算锂离子沿 <chem>LiFePO4</chem> 不同晶向的扩散速率。教程主要包括如下内容：

导入 <chem>LiFePO4</chem> 晶体结构
优化 <chem>LiFePO4</chem> 晶格参数
创建 Li$_{1-x}$FePO$_4$ 结构
优化初始、末态结构
创建 NEB 轨迹（IDPP插值）
优化 Li 扩散路径
使用简谐近似过渡态理论 (HTST) 计算反应速率

提示

本教程使用特定版本的QuantumATK创建，因此涉及的截图和脚本参数可能与您实际使用的版本略有区别，请在学习时务必注意。

导入 LiFePO4 晶体结构

VNL 内部数据库没有 <chem>LiFePO4</chem> 的晶体结构。您可以从其他数据库下载结构文件，一般而言 CIF 文件通用性比较强。

方法1. 直接从Crystallography Open Database检索、下载 CIF 文件，导入 VNL。

下载并拷贝 CIF文件到 Project 文件夹，这样 LabFloor 将会自动加载它。其它方式，例如导入XYZ、CIF、CAR、VASP文件，参考VNL指南。

方法2. 从2016版开始，VNL图形用户界面可以直接检索COD数据库。在 Builder 界面里，选择 Add → from plugin → Crystallography Open Database。

优化LiFePO4晶格参数

将 <chem>LiFePO4</chem> 导入到 Builder 中，可以进一步仔细查看晶体结构。

注意：

为了更好地与参考文献^[4]对比，可以使用Bulk Tools ‣ Swap Axes ‣ C↔A和Z↔X交换晶向和坐标轴，使其与文献一致。但这样的话，结构会失去对称性类型的信息，本例中会变成简单正交晶系。使用 Bulk Tools ‣ Lattice Paramters 将 Lattice type 设置为 Simple orthorhombic。

为了优化晶格参数，将 <chem>LiFePO4</chem> 结构从 Stash 送入到 Scripter 中设置新的 DFT 计算任务，如下：

1，添加New Calculator，并修改如下参数：

SGGA-PBE 交换相关势
7x5x3 k-point sampling

（我们将使用默认的FHI赝势以及DZP基组）。

2，添加OptimizeGeometry，并修改：

Max forces：0.02 eV/Å；
Max stress：0.001 eV/Å3；
去掉 Constrain cell的☑️，允许晶格的优化；
勾选 Save trajectory，并修改名字为 LiFePo4_opt.nc

3，将默认输出文件名改为 LiFePo4_opt.nc。

最后将脚本送入 JobManager 并运行该计算。

用 Builder 打开优化之后得到的名为 glD002 的 BulkConfiguration，检查优化之后的晶格参数与实验数值仅有 0.1% 的差别。

创建Li(1-x)FePO4结构

本文中将会研究 Li 缺陷 Li$_{1-x}$FePO$_4$ 中 Li 原子的扩散。具体而言，将会删除晶体结构中四个 Li 原子中的一个。应注意到，扩散能垒与Li 浓度有关^[2]。

在下一节，设置 Li 原子沿着两个不同方向扩散的计算，因此为同一个 Li 原子需要创建一个初始结构和两个末态结构分别沿着 B、C 方向扩撒。

创建沿 B 方向扩散的初、末态结构

1，右键点击 Stash 中的结构，并复制两次这样就有了三个一摸一样的结构。

2，从第一个结构中，像下图所示一样，删除一个 Li 原子，作为 NEB 计算的初始结构。可以将其重命名为 initial。

3，双击激活第二个结构，之后重命名为 final_1。选择与初始结构中删除的相同的 Li 原子，打开Periodic Table 工具，点击金元素，从而将 Li 原子替换为 Au。

小贴士：

这里是在用一个小技巧将 Li 原子移动到特定位置。在本例中，Li 空位到位置通过这个 Au 原子来确定了。实际上也有几种其它的方法可以创建末态结构。例如，通过拷贝删除掉的那个 Li 原子的坐标。非常欢迎大家在这里是用其它相关的，你喜欢的方式在此处使用。

4，选择 27 号 Li 原子（也就是近原子沿着 B 方向相邻的那个原子），使用 Move 工具，把它移动到 Au 原子上（重合），如下图所示，确保 Snap 被勾选了。

5，最后，使用 Fuse 选项，删除掉 Au 原子，这样末态结构就好了。

创建沿C方向扩散的初、末态结构

用同样的过程，创建第二个末态结构，如下图所示（本例中是 24 号 Li 原子）

优化初、末态结构

将每一个结构：initial、final_1、and final_2 送到 Scripter 中，三者使用与 <chem>LiFePO4</chem> 完全相同的参数设置结构优化计算。不过本例中不要优化晶格参数，因此勾选 Constrain cell，如下图所示：

另外，本例中可以使用 unpolarized GGA.PBE 泛函。这样计算速度会快很多，对本例而言，不影响结果的正确性。这里可以下载三个输入脚本:initial.zip、final_1.zip和final_2.zip。

创建NEB路径的初始猜测

NEB方法研究Pt在Pt表面的扩散教程中，演示了如何设置NEB，这里只将主要步骤列出来。

1，末态结构都优化完成之后，找到LabFloor加载进来的优化后的结构gID002，将两个BulkConfiguration直接拖到Builder里面。

2，在Builder中打开Builders ‣ Nudged Elastic Band，并从Stash中分别将初始结构和对应B方向扩散的末态结构，也就是initial.nc和final_1.nc分别拖入其中。

3，选择Image Dependent Pair Potential插值方法，使用包括末态结构在内的9个映像。

4，使用initial.nc和final_2.nc创建Li沿着C方向扩散的初始路径。

优化 Li 扩散路径

现在可以设置并优化 L i在 <chem>LiFePO4</chem> 中的扩散路径。

1，分别将两个 NEB任务从 Stash 送到 Scripter 中。

2，添加并设定一个 New Calculator，参数选择与末态结构计算保持一致。

3，添加一个 OptimizeGeometry，并确保选择了 Climbing Image Method，因为需要用这个方法去寻找正确的过渡态，用于简谐近似过渡态理论分析（下文有介绍）。

4，最后运行两个 NEB 计算 neb_b.zip 和 neb_c.zip。注意，如同 NEB方法研究Pt在Pt表面的扩散中描述的那样，可以使用 MPI 多进程并行计算。尤其是 NEB 计算各个image的时候，各个进程本身就是并发的。QuantumATK 会尝试检测最佳的并行化策略，并在 log 文件中输出一小段相关的报告。本例中能够看到程序检测到一个关于负载均衡的问题。可以修改 MPI 进程数，从而达到最佳性能。

+------------------------------------------------------------------------------+
|                                                                              |
| Atomistix ToolKit 2015.1 [Build 5b36b25]                                     |
|                                                                              |
+------------------------------------------------------------------------------+
+------------------------------------------------------------------------------+
| NEB image level parallelism enabled.                                         |
+------------------------------------------------------------------------------+
| Total number of interior images: 7                                           |
| Total number of processes: 16                                                |
| Processes per image: 2                                                       |
| Images calculated in parallel: 8                                             |
|                                                                              |
| Process group # 0 will calculate images: 1                                   |
| Process group # 1 will calculate images: 2                                   |
| Process group # 2 will calculate images: 3                                   |
| Process group # 3 will calculate images: 4                                   |
| Process group # 4 will calculate images: 5                                   |
| Process group # 5 will calculate images: 6                                   |
| Process group # 6 will calculate images: 7                                   |
| Process group # 7 will be idle!                                              |
| WARNING: Load balance problem detected.                                      |
|          The number of interior images calculated in parallel is not a       |
|          divisor of the number of interior images. This will greatly         |
|          reduce parallel performance. Either modify the parameter            |
|          "processes_per_image" or the number of MPI processes.               |
|                                                                              |
|          With 2 processes_per_image ideal performance will be obtained       |
|          by using 14 MPI processes.

小贴士

通过设置比较大的力收敛阈值，比如 0.08 eV/Å，这样计算会使用简单 NEB 计算，从而很快就会完成，不会启用 Climbing Image 方法。这样一般能够更快地给出一个能垒和扩散路径的粗略猜测。之后，可以设置更精确的 NEB 计算，并使用 Climbing Image 方法完成结果。

分析结果

只要计算正确收敛，就能在 LabFloor 中找到加载的内容中找到优化后的 NudgedElasticBand。

小贴士

计算结束后，一定不要忘记检查log文件，尤其是要检查是否收敛：

+------------------------------------------------------------------------------+
| NEB Optimization using the LBFGS optimizer                                   |
+------------------------------------------------------------------------------+
| Iteration   Step Length    Max. Force  Max. Energy Image    Max. Energy      |
+------------------------------------------------------------------------------+
|         0     0.000e+00     1.473e+00                  4       0.925917      |
|         1     4.825e-02     1.362e+00                  4       0.794101      |
...
|        23     1.220e-02     2.772e-01                  4       0.405770      |
|        24     3.731e-03     4.938e-02                  4       0.404796      |
+------------------------------------------------------------------------------+
| NEB optimization converged after 24 steps.                                   |
+------------------------------------------------------------------------------+

下面两个图中，使用 Movie Tool 显示了优化之后的 Li 原子沿着 B 和 C 方向扩散过程。

发现 Li 原子沿 C 方向移动需要跨越一个非常高的能垒，2.3eV；而沿着 B 方向，一个 Li 原子一个通道，移动则容易得多。这些结果与参考文献一致^[2]。

使用简谐近似过渡态理论计算反应速率

关于如何计算反应速率的细节，以及简谐近似过渡态理论本身，请参考 用简谐近似下的过渡态理论研究反应速率。

分析过程如下：

1，打开 Scripter。

2，双击 Analysis from file，选择 neb_B.nc 中的 glD001。

3，添加 HTSTEvent 分析，并设置 prefactor为 1e+13 1/s。

4，运行分析任务，一分钟左右完成。脚本文件：HTST_analysis.py

为了将 HTSTEvents 分析结果可视化，选择 LabFloor 中对应的对象，选中窗口右侧的 HTST Rates 插件。

这里，可以为扩散机制计算正向和你想反应速率常数。本例中，两个值非常接近，因为初始结构和末态结构基本上是等价的。

这里，也可以设置 Arrhenius 点图选项：

下表总结了本例中的结果：

Direction	Barrier	$k_{HTST}$
B	0.41	$1.6\times 10^{6}$
C	2.31	$1.5\times 10^{-26}$

参考文献

[1] M. Saiful Islam and Craig A. J. Fisher. Lithium and sodium battery cathode materials: computational insights into voltage, diffusion and nanostructural properties. Chem. Soc. Rev., 43(1):185, 2014. doi:10.1039/C3CS60199D.
[2] Chuying Ouyang, Siqi Shi, Zhaoxiang Wang, Xuejie Huang, and Liquan Chen. First-principles study of li ion diffusion in lifepo 4. Physical Review B, 2004. doi:10.1103/PhysRevB.69.104303.
[3] William Davidson Richards Yan Wang. Design principles for solid-state lithium superionic conductors. Nature Materials, 14(10):1026, 2015. doi:10.1038/nmat4369.
[4] Dhamodaran Santhanagopalan Yuri Janssen. Reciprocal salt flux growth of lifepo 4 single crystals with controlled defect concentrations. Chemistry of Materials, 25(22):4574, 2013. doi:10.1021/cm4027682.
英文教程原文：http://docs.quantumwise.com/tutorials/li_ion_diffusion.html

目录