致谢：
原notebook链接：https://bohrium.dp.tech/notebooks/73117386918
一、碘的基本性质
碘（化学符号：I，原子序数：53）是一种化学元素，在化学元素周期表中位于第5周期，第ⅦA族，属于卤素族。基态电子排布为。
碘于1811年由法国化学家贝尔塔·路易斯·吕瓦尔首次发现和分离。碘是一种深紫色的固体，在常温常压下呈片状晶体。它具有较强的挥发性，当受热时，固体碘会直接转变为紫色的气体，这个过程称为升华，在水中，碘不溶解，但可以溶解于一些有机溶剂中，例如乙醇和二甲苯。碘化学性质活泼且具有氧化性，在化合物中，碘主要以-1价存在，此时常见碘化物包括碘化氢、碘化钠、碘化钾、四碘化碳、碘化银和三碘化氮等，碘及其化合物在生物学、医学、工业和其他领域中有多种重要作用和用途。
二、DFT算法简介
密度泛函理论（DFT）是一种研究多电子体系电子结构的方法。密度泛函理论在物理和化学上都有广泛的应用，特别是用来研究分子和凝聚态的性质，是凝聚态物理计算材料学和计算化学领域最常用的方法之一。
密度泛函理论的概念起源于Thomas-Fermi模型，但直到Hohenberg-Kohn定理提出之后才有了坚实的理论依据。Hohenberg-Kohn第一定理指出体系的基态能量仅仅是电子密度的泛函。Hohenberg-Kohn第二定理证明了以基态密度为变量，将体系能量最小化之后就得到了基态能量。
密度泛函理论最普遍的应用是通过Kohn-Sham方法实现的。 在Kohn-Sham的框架中，最难处理的多体问题（由于处在一个外部静电势中的电子相互作用而产生的）被简化成了一个没有相互作用的电子在有效势场中运动的问题。这个有效势场包括了外部势场以及电子间库仑相互作用的影响，例如，交换相关作用。处理交换相关作用是KS方法中的难点。并没有精确求解交换关联势的方法。最简单的近似求解方法为局域密度近似(LDA近似)。LDA近似使用均匀电子气来计算体系的交换能（均匀电子气的交换能是可以精确求解的），而相关能部分则采用对自由电子气进行拟合的方法来处理。其他交换关联势包括广义梯度近似（GGA）等。
更多有关DFT计算原理的内容，参见了解第一性原理计算与密度泛函理论
三、ABACUS计算方法介绍
3.1 ABACUS电子自洽迭代计算（SCF）
ABACUS 是一款基于密度泛函理论（DFT）的第一性原理计算软件，自洽场计算（SCF, self-consistent field calculation）是其计算的核心过程。在某些教程中，这个过程也被称为「电子步」，即通过优化电荷密度来寻找势能最低点的过程，在 DFT 计算中，首要任务是分析体系的势能面。通过 SCF 计算，可以得到一个体系的基态结构和基态能量。
ABACUS支持LCAO（Linear Combination of Atomic Orbital，原子轨道线性组合）基组和PW（平面波）基组进行SCF计算，本教程中采取PW基组进行计算。 关于SCF计算的更多细节，参见ABACUS 使用教程｜电子自洽迭代
3.2 ABACUS结构优化（晶格弛豫）
在第一性原理计算中，晶格弛豫（Lattice relaxation）指的是一种过程，即通过调整晶格内原子的位置和晶格参数（如晶格常数），使晶体结构达到最低能量状态。
晶格弛豫的主要目的是找到晶体结构的稳定构型，这对于预测和理解材料的性能至关重要。在这个过程中，计算机模拟会根据所采用的理论框架来迭代地优化原子位置和晶格参数，最终得到一个能量最低的结构。这个结构通常被认为是晶体在实际条件下所处的状态。
具体的计算条件与SCF计算类似，更多参见ABACUS 使用教程｜结构优化
四、计算步骤和数据处理
4.1 晶体收敛性测试
由于ABACUS3.6.0的pw计算不支持uspp赝势，使用的是官网给出的SG15-v1.0赝势。 关于输入文件，查阅得到元素单质晶体为hcp结构，stru模板可以对照上文链接中算例写出，碘晶体a=4.934Å，c=6.695Å。而input和kpt参数改变不多，后面着重讨论参数改变时的收敛性。
首先是关于ecut的收敛性，平面波（pw）作为一种可以用于描述周期性边界条件下的电子波函数和电荷密度基矢量，它们都是正交的，且可以通过一个 ecutwfc 参数来控制基矢量的个数，ecut 其实代表了每一个平面波所对应的动能，如果 ecut 取得越大，则基矢量可以描述震荡得越剧烈的物理量（例如氧原子的 2p 轨道），那么计算结果就会越精确，但同时所带来的机时成本消耗也越大。因此我们采用 pw 基组做真正计算前，需要对 ecut 进行测试来获得一个足够准确且效率也高的取值。默认k值为4（一共444），将ecut值从20取到100Ry，得到晶体的体系总能量如下图所示：
可以看到，约在ecut=50Ry时，就可以认为体系总能量收敛了（之后变化小于0.01eV），收敛性很好。
然后讨论关于k的收敛性，取不同的k值（总节点数为k的立方），由上问结论，取ecut=50Ry已经足够，k从2取到8。得到计算时间t随k变化如下图所示：
体系总能量E随k的变化如下图所示：
可以看出，计算时间t随着k的增加会显著上升，而晶体能量E随着k的增加呈现出一种宏观上的收敛趋势，受制于测量的设置，这种趋势可能并不是十分严格，但是有理由认为在k趋于无穷时，总能量E能够达到收敛。针对ecut和k的收敛性测试，在一定程度上印证了该计算模型和方法（DFT）的可行性。
4.2 晶格结构优化计算
首先对上文所述的碘晶体hcp结构进行优化，参数与上文3.2中设置相同，ecut=100Ry，应力收敛阈值2kbar，k=4。设置初始值为a=4.934Å，c=6.695Å，空间群P6/mmm(191)进行晶格结构优化，迭代数次便收敛得到结果a=5.529Å，c=4.605Å。
实验没有查到相同空间群的晶格结构，有类似的晶体结构a=5.117Å，但是c=9.832Å，考虑到实验中这种晶体结构的γ=126.540°，而非严格的120度，应该将DFT计算得到的c乘2得到c=9.210Å，DFT计算值与实验值比较接近，可以认为结构优化结果较为合理。
然后以碘化钾（KI）晶体为例进行优化计算，碘化钾晶体结构类似氯化钠，取初始晶格常数a=4.436Å，空间群Fm-3m(225)进行晶格结构优化，也能比较快的得到收敛结果a=6.915Å，晶格矢量均正交。
实验数据为a=5.081Å,但是晶格夹角为60度，将计算得到a除以,得到取这种晶胞时的晶格常数a=4.890Å，与实验值相当接近，也可以认为是合理的。
最后分析误差来源，个人认为主要来自模型和计算误差，碘原子的价电子过多，赝势的描述可能不是特别精确，笔者的晶体建模与实际也有所出入，导致碘单质晶体计算结果不是很好。