CASTEP计算理论总结+实例分析

CASTEP 计算理论总结

XBAPRS

CASTEP 特点是适合于计算周期性结构，对于非周期性结构一般要将特定的部分作为周期性结构，建立单位晶胞后方可进行计算。CASTEP 计算步骤可以概括为三步：首先建立周期性的目标物质的晶体；其次对建立的结构进行优化，这包括体系电子能量的最小化和几何结构稳定化。最后是计算要求的性质，如电子密度分布(Electron density distribution)，能带结构(Band structure)、状态密度分布(Density of states)、声子能谱(Phonon spectrum)、声子状态密度分布(DOS of phonon)，轨道群分布(Orbital populations)以及光学性质(Optical properties)等。本文主要将就各个步骤中的计算原理进行阐述，并结合作者对计算实践经验，在文章最后给出了几个计算事例，以备参考。

CASTEP 计算总体上是基于DFT ，但实现运算具体理论有：

离子实与价电子之间相互作用采用赝势来表示；

超晶胞的周期性边界条件；

平面波基组描述体系电子波函数；

广泛采用快速fast Fourier transform (FFT) 对体系哈密顿量进行数值化计算；

体系电子自恰能量最小化采用迭带计算的方式；

采用最普遍使用的交换-相关泛函实现DFT 的计算，泛函含概了精确形式和屏蔽形式。 一， CASTEP 中周期性结构计算优点

与MS 中其他计算包不同，非周期性结构在CASTEP 中不能进行计算。将晶面或非周期性结构置于一个有限长度空间方盒中，按照周期性结构来处理，周期性空间方盒形状没有限制。之所以采用周期性结构原因在于：依据Bloch 定理，周期性结构中每个电子波函数可以表示为一个波函数与晶体周期部分乘积的形式。他们可以用以晶体倒易点阵矢量为波矢一系列分离平面波函数来展开。这样每个电子波函数就是平面波和，但最主要的是可以极大简化Kohn-Sham 方程。这样动能是对角化的，与各种势函数可以表示为相应Fourier 形式。

```2[()()()]``,,k G V G G V G G V G G C C ion H xc i i k G GG i k G δε∑++-+-+-=++ 采用周期性结构的另一个优点是可以方便计算出原子位移引起的整体能量的变化，在CASTEP 中引入外力

或压强进行计算是很方便的，可以有效实施几何结构优化和分子动力学的模拟。平面波基组可以直接达到有效的收敛。

计算采用超晶胞结构的一个缺点是对于某些有单点限缺陷结构建立模型时，体系中的单个缺陷将以无限缺陷阵列形式出现，因此在建立人为缺陷时，它们之间的相互距离应该足够的远，避免缺陷之间相互作用影响计算结果。在计算表面结构时，切片模型应当足够的薄，减小切片间的人为相互作用。 CASTEP 中采用的交换-相关泛函有局域密度近似（LDA ）（LDA ）、广义梯度近似（GGA ）和非定域交换-相关泛函。CASTEP 中提供的唯一定域泛函是CA-PZ ，Perdew and Zunger 将Ceperley and Alder 数值化结果进行了参数拟和。交换-相关泛函的定域表示形式是目前较为准确的一种描述。

采用梯度校正的非定域或广义梯度近似泛函与电子密度梯度d dr ρ和电子密度ρ都有关，这样可以同时提高能量和结构预测的准确性，但计算耗时。CASTEP 中提供的非定域泛函有三种：PBE 泛函与PW91泛函计算在本质上实际是相同的，但在电子密度变化迅速体系中PBE 泛函实用性更好；RPBE 是特别用来提高DFT 描述金属表面吸附分子能量的泛函，White and Bird 描述了各种梯度校正泛函计算方法，利用广义梯度近似计算总能量使用平面波基组与定域泛函相比并不直接。包含梯度近似的交换-相关泛函计算时对电子密度数据的精度要求较高，对计算机内存占用会增大。通过采用与平面波基组总能量计算中分裂交换-相关能量采用一系列空间网格相一致的方法来定义交换-相关势。

平面波基组（Plane wave basis set ）

Bloch 理论表明每个k 点处电子波函数都可以展开成离散的平面波基组形式，理论上讲这种展开形式包含的平面波数量是无限多的。然而相对于动能较大的情况，动能|k +G |2

很小时平面波系数C k +G 更重要。Name Description Reference

PW91 Perdew-Wang generalized-gradient approximation, PW91 Perdew and Wang

PBE Perdew-Burke-Ernzerhof functional, PBE Perdew et al.

RPBE Revised Perdew-Burke-Ernzerhof functional, RPBE

Hammer et al.

调节平面波基组，其中包含的平面波动能小于某个设定的截止能量，如图所示（球体半径与截止能量平方根成比例）：

总能量计算会因为平面波特定能量截止而产生误差，通过增加体系能量截止数值就可以减小误差幅度。理论上截止能量必须提高到总能量计算结果达到设定的精确度为止，如果你在进行关于相稳定性的研究，而需要对比每个相能量的绝对值时，这是一种推荐计算方法。不过，同一个结构在低的截止能量下收敛引起的差别要小于总体能量本身。因此可以选用合适的平面波基组对几何结构进行优化或进行分子动力学研究。以上的方法对Brillouin区取样收敛测试同样成立。

有限平面波基组的校正

采用平面波基组的一个问题是截止能量与基组数量的变化是间断的，一般而言在ｋ点基组（k-point set）中不同ｋ点对应不同能量截止（cutoffs）时就会产生这种不连续性。此外，在截止能量不变时，晶胞形状和尺寸的变化都会引起平面波基组的间断。通过采用更加致密的ｋ点基组就可以解决这个问题，与特定平面波基组相关的加权性也会消除。然而即使在ｋ基组取样很致密的情况下，这个问题依然存在，对其近似的解决方法就是引入一个校正因子（correction factor），利用某个状态基组计算使用了无限数量的ｋ点与实际采用的数量之间的差别来确定。晶体结构在进行几何优化时如果基组不能真正的达到绝对的收敛，有限基组的纠正就很重要。比如硅的规范－保守赝势很“软”，在平面波基组截止能量是２００ｅＶ时就已经可以得到准确的计算结果了。但如果计算状态方程时使用上述截止能量（比如体积与总能量和压强都有关系），能量最小时对应的体积与体系内压为零时对应体积是不同的。在提高截止能量和增加ｋ点取样基础上重复对状态方程的计算，这两个体积之间的差别会越来越小。此外截止能量低时计算得到的E-V曲线呈现锯齿状，提高截止能量计算的曲线连续而平滑。E-V曲线中出现锯齿状的原因在于平面波基组在相同的截止能量时由于晶体点阵常数不同引起的平面波基组数量的间断。对总能量进行有限基组的校正，使得我们可以在一个恒定数量基组状态下进行计算，即使采用了恒定的截止能量这个更强制条件也可以纠正计算结果。Milman等详细的讨论了这种计算方法的细节。进行这种校正所需要的唯一的参数就是dE tot/d lnE cut，E tot是体系总能量，E cut是截止能量。dE tot/d lnE cut的值很好的表示了能量截止和ｋ点取样计算收敛性质。当它的数值（每个原子）小于 eV/atom时，计算就达到了良好的收敛精度，对于大多数计算 eV/atom就足够了。