这是本文档旧的修订版！
Li-S 电池的开路电压：ReaxFF 方法分子动力学
版本：2017.1
硫是用于可再充电能量存储设备中较为有前途的阴极材料。电池放电期间，锂离子从富锂阳极迁移到硫阴极，并与硫反应（阴极锂化）：
$$x \textrm{Li} + \textrm{S} \rightarrow \textrm{Li}_x \textrm{S}$$
开路电压（OCV）是该锂化过程中电化学驱动力。
在本教程中，我们将利用分子动力学（MD）模拟来计算电池放电期间的 OCV 曲线。我们采用了 ATK-ForceField 计算器和在参考文献 [IOB+15] 中提到并由伊斯兰教同事开发的 ReaxFF 力场。模拟退火法可用于模拟一系列成分为 x 的无定形 LixS 化合物，并将所得 OCV 曲线与参考文献 [IOB+15] 中的图 3 进行比较。
注意
参考文献 [IOB+15] 中的图 3 是由采用混合了大规模典型的蒙特卡洛/分子动力学（GC-MC / MD）方法获得的数据制作而成，该方法目前在 QuantumATK 中尚无法实现。我们没有使用 GC-MC 将 Li 原子置于低能位，而是在这里使用了一种模拟退火的方法。
开路电压的理论定义
忽略焓和熵对能量的贡献，OCV 根据总能量差计算得来：
$$\textrm{OCV} = -\frac{E_{\textrm{Li}_x \textrm{S}} - nE_\textrm{Li} - mE_\textrm{S}}{n} ,$$
在这里，$n$($m$) 是 LixS 中锂（硫）原子的数量（如$x=n/m$），$E_\textrm{Li}$ 和 $E_\textrm{S}$ 分别为每个纯锂和硫晶体的总能量。
概要
我们首先考虑 Li0.4S，并描述为计算纯 S 部分锂化到 Li0.4S 的开路电压的所需步骤：
- 无定形 Li0.4S 化合物
- 模拟退火
- 开路电压
然后利用 ATK 的 Python 脚本自动模拟更多范围的 LixS 化合物。然后计算得到完整的 OCV 曲线，并研究不同成分的径向分布函数：
- 完整的开路电压图
- 径向分布函数
提示
本教程提供了所有可供下载的 QuantumATK Python 脚本。您也可以在此处下载压缩文件：↓ scripts.zip。
注意
这是一个高级教程，因此不会详细地介绍如何操作 ATK-ForceField 计算器和 QuantumATK 分析工具的相关基础知识。如有需要，请参阅教程 Molecular Dynamics: Basics。
无定形 Li0.4S 化合物
在电池工作条件下，阴极材料不会如期待的保持完美的晶相；它最可能是无定形的。QuantumATK 提供了多种用于生成非晶结构的工具，例如 Amorhpous Prebuilder 和 PackMol 工具。本教程中我们使用后者。
打开 QuantumATK，创建一个新的项目，最好放在硬盘驱动的新目录中。然后打开 
Builder。
PackMol
我们使用 PackMol 生成无定形 Li0.4S 化合物的初始结构，它可以将输入分子构型随机混合进一个合理的无定形块体构型。因此，您需要先向 Builder 的 Stash 添加硫原子和锂原子：
- 在 Stash，点击 Add
New Configuration，在分子构型中添加一个氢原子。 - 选中氢原子，利用
工具将 H 原子转变为 S。 - 重复上述两个步骤，创建锂原子。
- 分别将 Stash 中的两项重命名为 “sulfur” 和 “lithium”。
转至 Builders 
PackMol 插件并点击，然后将打开 PackMol 程序。然后按照以下步骤创建包含 2048 个硫原子的 Li0.4S 无定形晶体。
- 将 Stash 中的硫原子和锂原子拖放到 PackMol 程序。构型将会被添加到 “Molecule type” 的列表。
- 将硫原子的数目增加到 2048，锂原子的数目输入为 819（由2048 × 0.4 = 819）。
- 其他设置保持默认，点击 Create。
无定形 Li0.4S 结构已被自动添加到 Stash。重命名为 “Li04S”。
结构优化
PackMol 产生的无定形结构可确保在所有混合实体之间具有最小距离（默认为 2 Å），但该结构可能仍远未达到平衡。因此，在使用 MD 将系统加热到高温之前，最好对原*子坐标和晶格矢量进行优化。
- 将 “Li04S” 构型发送到
Script Generator，添加 
New Calculator 和 
OptimizeGeometry 脚本模块。 - 打开 calculator 模块，选择 ATK-ForceField 计算器。自动选择的力场为 ReaxFF_LiS_2005 —— 在参考 [IOB+15] 中生成。
- 打开 geometry optimization 模块，设置 force tolerance 和最大优化步数分别为 0.5 eV/Å 和 1000。然后勾选 Save trajectory 选项，输入所保存弛豫轨迹文件名称为
Li04S。
保存 QuantumATK Python 脚本为 Li04S.py 并使用 
Job Manager 或在终端执行以下命令运行：
atkpython Li04S.py 2>&1 | tee Li04S.log
您也可以在此处下载脚本：↓ Li04S.py。
重要
相较于 MPI，使用 OpenMP 线程可以最有效地在多个物理内核上进行 ATK-ForceField 并行计算，请参阅教程中关于使用 Job Manager 仅通过线程进行https://docs.quantumatk.com/tutorials/job_manager_local/job_manager_local.html#job-manager-local-mkl本地和远程 QuantumATK 计算的部分。
如果从终端运行 QuantumATK，则应取消设置环境变量 OMP_NUM_THREADS 或将其设置为您打算使用的最大 OpenMP 线程数。我们建议将 MKL_DYNAMIC = TRUE 设置为允许 MKL 动态地更改线程数。
结构优化现在已经结束了，输出数据应在保存在 QuantumATK 的 LabFloor：
您可以利用 
Viewer 将弛豫轨迹可视化，您会看到，在最小化原子力和应力的作用下，晶胞膨胀了一点。
模拟退火
无定形 Li0.4S 化合物的室温结构（原子坐标和晶胞矢量）不是事先知道的，因此我们使用模拟退火的方法找到它。基本思想是将系统加热到 1600 K，然后使用分子动力学将其缓慢冷却至 300 K，所有这些过程都在固定压力下进行，但允许晶胞体积调整（NPT 动力学）。
计算流程：
- 让系统在 1600 K 达到平衡。
- 让系统缓慢地从 1600 K 冷却至 300 K。
- 让系统在 300 K 达到平衡。
在 1600K 达到平衡
在 LabFloor （保存在 Li04S.hdf5）标签为 BulkConfiguration_1 的块体结构即为优化后的。将其选中并拖放到 
Scripter，会自动地添加适合的 ATK-ForceField 计算器模块，因此添加的是 
MolecularDynamics 脚本模块。
首先，设置默认输出文件名称为 Li04S_eq-1600K.hdf5：
打开 MD 脚本模块，调整一下设置：
- 设置 MD 类型为 NPT Berendsen。
- 增加 MD 步数至 50000。
- 增加 log interval 为 1000。
- 将 Reservoir temperature 和 Final temperature 都设置为 1600 K。
- 将初始粒子速度 Maxwell–Boltzmann 分布的温度也设置为 1600 K。
- 选项：不勾选 Save trajectory 选项。
- 点击 OK，保存设置。
保存脚本为 Li04S_eq-1600K.py 并运行。取决于计算资源（即 OpenMP 线程数），该作业可能需要一个小时才能完成。
您也可以在此处下载 QuantumATK Python 脚本：↓ Li04S_eq-1600K.py。
一旦计算完成，最终的 MD 轨迹应在 LabFloor 上弹出。可以使用 Movie Tool 将 MD 模拟进程可视化。如果需要，您可以使用鼠标右键单击启用 Show bonds 选项。请注意，在整个模拟过程中温度和能量如何大致保持恒定，且晶胞在开始时会明显膨胀，但随后原子在内部移动时保持相当稳定。因此构型在恒定温度下保持平衡。
