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

CASTEP 计算理论总结XBAPRSCASTEP 特点是适合于计算周期性结构，对于非周期性结构一般要将特定的部分作为周期性结构，建立单位晶胞后方可进行计算。

CASTEP 计算步骤可以概括为三步：首先建立周期性的目标物质的晶体；其次对建立的结构进行优化，这包括体系电子能量的最小化和几何结构稳定化。

最后是计算要求的性质，如电子密度分布(Electron density distribution)，能带结构(Band structure)、状态密度分布(Densityof states)、声子能谱(Phonon spectrum)、声子状态密度分布(DOS of phonon)，轨道群分布(Orbitalpopulations)以及光学性质(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 数值化结果进行了参数拟和。

交换-相关泛函的定域表示形式是目前较为准确的一种描述。

Name DescriptionReference PW91 Perdew-Wang generalized-gradient approximation, PW91Perdew and Wang PBE Perdew-Burke-Ernzerhof functional, PBEPerdew et al. RPBE Revised Perdew-Burke-Ernzerhof functional, RPBE Hammer et al.采用梯度校正的非定域或广义梯度近似泛函与电子密度梯度ddrρ和电子密度ρ都有关，这样可以同时提高能量和结构预测的准确性，但计算耗时。

CASTEP中提供的非定域泛函有三种：PBE泛函与PW91泛函计算在本质上实际是相同的，但在电子密度变化迅速体系中PBE泛函实用性更好；RPBE是特别用来提高DFT描述金属表面吸附分子能量的泛函，White and Bird描述了各种梯度校正泛函计算方法，利用广义梯度近似计算总能量使用平面波基组与定域泛函相比并不直接。

包含梯度近似的交换-相关泛函计算时对电子密度数据的精度要求较高，对计算机内存占用会增大。

通过采用与平面波基组总能量计算中分裂交换-相关能量采用一系列空间网格相一致的方法来定义交换-相关势。

平面波基组（Plane wave basis set）Bloch理论表明每个k点处电子波函数都可以展开成离散的平面波基组形式，理论上讲这种展开形式包含的平面波数量是无限多的。

然而相对于动能较大的情况，动能|k+G|2很小时平面波系数C k+G更重要。

调节平面波基组，其中包含的平面波动能小于某个设定的截止能量，如图所示（球体半径与截止能量平方根成比例）：总能量计算会因为平面波特定能量截止而产生误差，通过增加体系能量截止数值就可以减小误差幅度。

理论上截止能量必须提高到总能量计算结果达到设定的精确度为止，如果你在进行关于相稳定性的研究，而需要对比每个相能量的绝对值时，这是一种推荐计算方法。

不过，同一个结构在低的截止能量下收敛引起的差别要小于总体能量本身。

因此可以选用合适的平面波基组对几何结构进行优化或进行分子动力学研究。

以上的方法对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的值很好的表示了能量截止和ｋ点取样计算收敛性质。

当它的数值（每个原子）小于0.01 eV/atom时，计算就达到了良好的收敛精度，对于大多数计算0.1 eV/atom就足够了。

非定域交换-相关泛函基于LDA或GGA的泛函的Kohn-Sham方程在计算能带带隙上存在低估。

这对晶体或分子相关性质以及能量的描述是没有影响的。

然而要理解半导体和绝缘体性质，就必须得到关于电子能带结构的准确的描述。

DFT能带带隙计算误差可以通过引入经验“剪刀”校正，相对于价带而言导带产生了一个刚性的变化。

当实验提供的能带带隙准确时，光学性质计算得到了较为准确的结果。

电子结构实验数据缺乏时采用“剪刀”工具进行预测性研究或对能带带隙调整是不可靠的。

关于DFT计算中能带带隙问题已经发展许多技术，但这些技术大多复杂而且很耗时，实际计算中最常用的是屏蔽交换（Sx-LDA），建立在广义Kohn-Sham方法基础上。

广义Kohn-Sham泛函允许我们将总能量交换分布泛函分离为非定域、定域以及屏蔽密度组元。

在CASTEP计算中采用的广义Kohn-Sham方法有：•HF: exact exchange, no correlation•HF-LDA: exact exchange, LDA correlation•sX: screened exchange, no correlation•sX-LDA: screened exchange, LDA correlation与LDA和GGA相比No local functionals 也有一些缺陷。

在屏蔽交换泛函中不存在已知形式应力张量表达方式，因此没有完全的非定域势可以用于单位晶胞结构优化或进行NPT/NPH动力学。

这样利用这些泛函计算的光学性质很有可能是不准确的。

在哈密顿量中引入一个完全非定域组元就可以解决这个问题，这个额外的矩阵元破坏了光学矩阵元素由位置算符转换为动量算符常用表达形式，使得哈密顿量对易很复杂。

规范保守赝势和超软赝势赝势是利用平面波基组计算体系总能量中关键的一个概念，价电子与离子实之间强烈的库仑势用全势表示时由于力的长程作用很难准确的用少量的Fourier变换组元表示。

解决这个问题的另一种方法从体系电子的波函数入手，我们将固体看作价电子和离子实的集合体。

离子实部分由原子核和紧密结合的芯电子组成。

价电子波函数与离子实波函数满足正交化条件，全电子DFT理论处理价电子和芯电子时采取等同对待，而在赝势中离子芯电子是被冻结的，因此采用赝势计算固体或分子性质时认为芯电子是不参与化学成键的，在体系结构进行调整时也不涉及到离子的芯电子。

为了满足正交化条件全电子波函数中的价电子波函数在芯区剧烈的振荡，这样的波函数很难采用一个合适的波矢来表达。

在赝势近似中芯电子和强烈库仑势替代为一个较弱的赝势作用于一系列赝波函数。

赝势可以用少量的Fourier变换系数来表示。

理想的赝势在芯电子区域是没有驻点的，因此需要平面波矢数量很少。

众所周知的是现在将赝势与平面波矢相结合对描述化学键是很有用的。

全离子势的散射性质可以通过构筑赝势得到重现，价电子波函数相位变化与芯电子角动量成分有关，因此赝势的散射性质就与轨道角动量是相关的。

赝势最普遍表达方式是：V= |lm> V l <lm|NLwhere |lm> are the spherical harmonics and V l is the Pseudopotential for angular momentum l.在不同角动量通道均采用同一个赝势值称为定域赝势（Local Pseudopotential），定域赝势计算效率更高，一些元素采用定域赝势就可以达到准确描述。