CASTEP计算原理---------XBAPRS
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形式。
?[k?G?2?V(G?G)?V(G?G)?V(G?G)]C??C``ionHxcii,k?GGGi,k?G```
采用周期性结构的另一个优点是可以方便计算出原子位移引起的整体能量的变化，在CASTEP中引入外力
或压强进行计算是很方便的，可以有效实施几何结构优化和分子动力学的模拟。平面波基组可以直接达到有效的收敛。
计算采用超晶胞结构的一个缺点是对于某些有单点限缺陷结构建立模型时，体系中的单个缺陷将以无限缺陷阵列形式出现，因此在建立人为缺陷时，它们之间的相互距离应该足够的远，避免缺陷之间相互作用影响计算结果。在计算表面结构时，切片模型应当足够的薄，减小切片间的人为相互作用。 CASTEP中采用的交换-相关泛函有局域密度近似（LDA）（LDA）、广义梯度近似（GGA）和非定域交换-相关泛函。CASTEP中提供的唯一定域泛函是CA-PZ，Perdew and Zunger将Ceperley and Alder数值化结果进行了参数拟和。交换-相关泛函的定域表示形式是目前较为准确的一种描述。
Name Description
PW91 Perdew-Wang generalized-gradient approximation, PW91 PBE Perdew-Burke-Ernzerhof functional, PBE
RPBE Revised Perdew-Burke-Ernzerhof functional, RPBE
Reference Perdew and Wang Perdew et al. Hammer et al.
1
CASTEP计算原理---------XBAPRS
采用梯度校正的非定域或广义梯度近似泛函与电子密度梯度
d?dr和电子密度?都有关，这样可以
同时提高能量和结构预测的准确性，但计算耗时。CASTEP中提供的非定域泛函有三种：PBE泛函与PW91泛函计算在本质上实际是相同的，但在电子密度变化迅速体系中PBE泛函实用性更好；RPBE是特别用来提高DFT描述金属表面吸附分子能量的泛函，White and Bird描述了各种梯度校正泛函计算方法，利用广义梯度近似计算总能量使用平面波基组与定域泛函相比并不直接。包含梯度近似的交换-相关泛函计算时对电子密度数据的精度要求较高，对计算机内存占用会增大。通过采用与平面波基组总能量计算中分裂交换-相关能量采用一系列空间网格相一致的方法来定义交换-相关势。
平面波基组（Plane wave basis set）
Bloch理论表明每个k点处电子波函数都可以展开成离散的平面波基组形式，理论上讲这种展开形式
2
包含的平面波数量是无限多的。然而相对于动能较大的情况，动能|k+G|很小时平面波系数Ck+G更重要。调节平面波基组，其中包含的平面波动能小于某个设定的截止能量，如图所示（球体半径与截止能量平方根成比例）：
总能量计算会因为平面波特定能量截止而产生误差，通过增加体系能量截止数值就可以减小误差幅度。理论上截止能量必须提高到总能量计算结果达到设定的精确度为止，如果你在进行关于相稳定性的研究，而需要对比每个相能量的绝对值时，这是一种推荐计算方法。不过，同一个结构在低的截止能量下收敛引起的差别要小于总体能量本身。因此可以选用合适的平面波基组对几何结构进行优化或进行分子动力学研究。以上的方法对Brillouin区取样收敛测试同样成立。
有限平面波基组的校正
采用平面波基组的一个问题是截止能量与基组数量的变化是间断的，一般而言在ｋ点基组（k-point set）中不同ｋ点对应不同能量截止（cutoffs）时就会产生这种不连续性。此外，在截止能量不变时，晶胞形状和尺寸的变化都会引起平面波基组的间断。通过采用更加致密的ｋ点基组就可以解决这个问题，与特定平面波基组相关的加权性也会消除。然而即使在ｋ基组取样很致密的情况下，这个问题依然存在，对其近似的解决方法就是引入一个校正因子（correction factor），利用某个状态基组计算使用了无限数量的ｋ点与实际采用的数量之间的差别来确定。晶体结构在进行几何优化时如果基组不能真正的达到绝对的收敛，有限基组的纠正就很重要。比如硅的规范－保守赝势很“软”，在平面波基组截止能量是２００ｅＶ时就已经可以得到准确的计算结果了。但如果计算状态方程时使用上述截止能量（比如体积与总能量和压强都有关系），能量最小时对应的体积与体系内压为零时对应体积是不同的。在提高截止能量和增加ｋ点取样基础上重复对状态方程的计算，这两个体积之间的差别会越来越小。此外截止能量
2
CASTEP计算原理---------XBAPRS
低时计算得到的E-V曲线呈现锯齿状，提高截止能量计算的曲线连续而平滑。E-V曲线中出现锯齿状的原因在于平面波基组在相同的截止能量时由于晶体点阵常数不同引起的平面波基组数量的间断。对总能量进行有限基组的校正，使得我们可以在一个恒定数量基组状态下进行计算，即使采用了恒定的截止能量这个更强制条件也可以纠正计算结果。Milman等详细的讨论了这种计算方法的细节。进行这种校正所需要的唯一的参数就是dEtot/d lnEcut，Etot是体系总能量，Ecut是截止能量。dEtot/d lnEcut的值很好的表示了能量截止和ｋ点取样计算收敛性质。当它的数值（每个原子）小于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变换系数来表示。理想的赝势在芯电子区域是没有驻点的，因此需要平面波矢数量很少。众所周知的是现在将赝势与平面波矢相结合对描述化学键是很有用的。全离子势的散射性质可以通过构筑赝势得到重现，价电子波函数相位变化与芯电子角动量成分有关，因此赝势的散射性质就与轨道角动量是相关的。赝势最普遍表达方式是：
3
CASTEP计算原理---------XBAPRS
VNL = ? |lm> Vl where |lm> are the spherical harmonics and Vl is the Pseudopotential for angular momentum l.在不同角动量通道均采用同一个赝势值称为定域赝势（Local Pseudopotential），定域赝势计算效率更高，一些元素采用定域赝势就可以达到准确描述。赝势的硬度（hardness）在赝势的应用中是一个重要的概念，当一个赝势可以用很少的Fourier变换组元就可以准确描述时称为“软赝势”，硬赝势与此相反。早期发展的准确规范保守赝势很快就发现在过渡元素和第一周期元素（C、N、O，等）中的描述十分“硬”，提高规范保守赝势收敛性质的各种方法都已经被提出，在CASTEP中采用了由Lin等提出的动能优化而来的规范保守赝势。 Vanderbilt提出了另一种更基本的方法，放宽规范保守赝势的要求，从而生成更软的赝势。在超软赝势方法中，芯电子区的赝平面波函数可以尽可能的“软”，这样截止能量就可以大幅度的减少。超软赝势与规范保守赝势相比除了“更软”以外还有其它的优点，在一系列预先设定的能量范围内遗传算法确保了良好的散射性质，从而使赝势获得更好变换性和准确性。超软赝势通常将外部芯区按照价层处理，每个角动量通道中的占据态都包含了复合矢。这样就增加了赝势的变换性和准确性，但同时是以消耗计算效率为代价的。可转移性是赝势的主要优点。赝势是通过孤立的原子或离子特定的电子排部状态下构建的，因此可以准确的描述原子在那些特定排部下芯区的散射性质。在相应条件下产生的赝势可以用于各种原子电子排部状态以及各种各样的固体中，同样也确保了在不同的能量范围内具有正确的散射状态。Milman给出了不同化学环境和一系列结构中采用赝势描述准确性事例。非定域赝势即使在最有效离散表示情况下，体系能量赝势计算依然占用了大量计算时间。此外，在倒易点阵空间采用非定域赝势会因原子数目增多而耗时以原子立方数增大，因此对于大体系是很适用的。赝势非定域性是指只有在超过原子芯区时它才会扩展，由于芯区是很小的，特别是当体系包含有许多的真空腔体时，在实空间采用赝势来计算就有很大的优势。这时计算量随体系中原子数目平方增长，因此是很适合大体系计算的。将电子划分为芯电子和价电子在处理交换-相关相互作用时会产生新问题，在原子芯区两个亚体系叠加在赝势产生过程中很难完全去屏蔽。在赝势能量算符中与电子密度存在非线形关系的项就是交换-相关能。Louie等采用了一种简单的方法来处理芯电子和价电子密度之间非线性的交换-相关能。这种方法在很大程度上提高了赝势的可变换性，特别是自旋极化的计算更为准确。当准芯区电子不能简单处理为价电子时非线性核校正就很重要。另一方面将他们简单地包含在价层亚体系中从本质上可以避免NLCC处理的必要性。 规范保守赝势： 采用赝势计算关键在于可以有效的对化学键的价电子进行可再现的近似，赝势与全势在超过离子实半径以后具有完全相同的函数形式。 Figure 1. Schematic representation of the all-electron and pseudized wave functions and potentials 两个函数平方幅度的积分数值应该是相同的，这等同于要求赝势波函数具有规范-保守性，比如每个赝波函数只能描述一个电子的行为。这样的条件就确保了赝势可以再现正确的散射 (Scattering Properties)性质。 生成赝势的典型方法如下所述：选择某个特定的电子排部状态（不一定就 4 CASTEP计算原理---------XBAPRS 是基态）全部电子计算在一个孤立的原子中进行。从而得到原子价电子能量本征值和价电子波函数。选择一个离子赝势或赝波函数参数形式，通过对参数的调节，使得赝原子计算和全电子原子赝势计算采用相同的交换-相关势，在超过截止半径后与价电子波函数形式相同，赝势的本征值等于价电子的本征值。如果电子波函数和赝势波函数满足正交归一，两者在截止半径以外的匹配性决定了规范-保守条件自动成立。离子赝势的截止半径是实际物理芯区的二到三倍。截止半径越小，赝势越“硬”而适用性（transferability）好。计算精度和效率决定了实际中采用的截止半径的大小。 规范-保守赝势优化 在固体计算中依据能量的截止存在一系列优化赝势的方法，Lin基于Rappe早期工作提出了下列赝势产生方法：在截止半径（cutoff radius）内，赝势波函数可以表达为： ?PSl4`(qR)?`(R)j(r)??aijl(qir),lic?lcjc(qiRc)?l(Rc)i?1 j(qr)li是球形Bessel函数，在r=0和r=Rc之间有（i-1）个零点。为保证赝势的实用性，截止半径越大越 好。超过截止矢量qc对动能最小化可得到系数?i。 ?qc2PS,?2PS2PS?E(qi,q)???dr?l(r)??l(r)??dqq?l(q)kc00?? 在第一个方程中让qc等于q4。其他的三个限制条件使得赝波函数在进行Lagrange连乘（Lagrange multipliers）时保持正交化（normalization），并且使赝波函数在Rc处的第一个二介偏微分是连续的。半径相关Kohn-Sham 方程反转标准步骤产生的一个具有理想收敛性质的平滑赝势函数。Lee提出了进一步改进的方法，在CASTEP数据库中固体规范保守赝势就是采用他的思想设计的。这种通用的方法消除了在特定的截止半径处赝波函数的二介偏微分必须是连续的条件，因为它是自动满足这个条件的。这样对于特定截止半径Rc允许我们通过调节qc提高赝势的精度和计算效率。 超软赝势（ULTRASOFT PSEDUPOTENTIAL） 为了能够使平面波基组计算中所采用的截止能量尽可能的小，Vanderbilt提出了超软赝势方法。众所周知规范-保守赝势在收敛优化中存在本身缺陷，所以就设计了另一种方法。超软赝势基础是在大多数情况下只有当紧密结合原子价轨道加权性分数大部分在芯区时，利用平面波基组计算才要求较高的截止能量。在这种情况下，减少平面波基组的唯一方法就是解除（violate）规范-保守赝势成立条件，将这些轨道中的电子从芯区移去。芯区的赝势就可以尽可能的“软”，从而使截止能量降低达到要求。从技术上讲，通过引入一个广义的正交归一化条件就可以完成。为了覆盖全部电子电荷，在芯区对由电子波函数模平方产生的电子密度进行适度放大（augmented）。电子密度划分成两部分：扩展在整个晶体中“软”部分和定域在芯区的“硬”部分。 5 CASTEP计算原理---------XBAPRS 固体中超软赝势公式 超软赝势中总能量与采用其他赝势平面波方法时相同，非定域势VNL表达如下： (0)IV??Dnm?nNlnm,I投影算符β和系数D 度可以表示为： （0） I?m 分别表征赝势和原子种类的差别，指数 I对应于一个原子位置 。总能量用电子密 ?2I(r)??In(r)????i(r)??Qnmini?nm,I??I???mi?? Φ是波函数， Q(r) 是严格位于芯区的附加函数(Augment function) 。超软赝势完全由定域部分， ion(0) Vloc(r) 和系数D, Q and ?确定，这些变量计算方法在下文中将做介绍。 引入一个广义正交归一条件来解除规范-保守赝势的限制条件： ?iS?j??ij S是哈密顿重叠算符（Hermitian overlap operator） I?S?1?qnm?nnm,II?m 系数q是通过对Q（r）积分得到，超软赝势的Kohn-Sham方程可以写为： H?iH代表了动能和定域势能之和，如下所示： ??S?ii 6 CASTEP计算原理---------XBAPRS I?IH?T?V??Dnmneffnm,Iion I?m 在Veff中包含离子定域势Vloc(r)，Hartree势和交换-相关势等项。通过定义一些新参数就可以将因附加 （augmented）电子密度而产生所有项全部包含在赝势的非定域部分。 DInm?D??drVeff(r)Q(r)nmnm(0)I 与规范-保守赝势对比，不同之处在于在超软赝势中存在重叠算符S，波函数与D有关而且事实上投影算 符函数β（projector function）数量要比规范-保守赝势中大两倍多。与附加（augmented）电荷相关的一系列计算可以在实空间(real space)中进行，这与函数中定域势的性质有关。多余的步骤不会对计算效率产生较大的影响。在Laasonen文献中提供了超软赝势计算的详细方法以及总能量微分表达式。 赝势生成：与规范-保守赝势情况一样，在自由原子上对所有的电子进行计算，得到屏蔽原子势VAE(r)（screened atomic potential）。每个角动量选择一系列的参考能量?l，一般两个能量参考点就足够了。这些能量参考范围必须包含良好散射性质，在每个参考能量处求解与半径相关的Kohn-Sham方程，得到规则初始点。选择截止半径，对上面产生的每个全电子波函数构筑一个赝势?，唯一的限制条件是它必须在Rcl处与?平滑相交。定义一个比所有芯区半径稍大的辅助半径R。最后就形成了定域轨道（超过R时就消失）： ?n?(??T?V)?nnlocB? 以及它们矩阵内积（inner products）：nmion ?n?m (0) 这样就可以定义用于固体计算的变量 (Vloc(r), D, Q and ?): Q(r)??(r)?(r)??(r)?(r);q??drQnm(r)nmnmnmnm?? 7 CASTEP计算原理---------XBAPRS ?n??(B?1)nm?mm ion (0) 采用去屏蔽（descreening procedure）方法计算Vloc(r), D系数: Vlocion(r)?V(r)?V(r)?V(r)locHxc (0)Dnm?Bnm??mqnm??drVloc(r)n(r)在去屏蔽方法中可以引入非线性核校正方法（The nonlinear core correction (NLCC)），这与规范-保守赝势中所采用的方法完全一致。在以下情况下超软赝势是很适用的： 赝本征值与所有电子本征值相同，在芯半径截止区以外赝轨道波函数与价电子波函数匹配一致；对于每个参考能量散射性质都是正确的，这样通过增加参考能量点数目就可以系统的提高赝势的适用性；在参考电子排部情况下，赝势价电子密度与全价电子密度相同。 关于非线性核校正 Louie等人第一次提出了非线性芯校正，使得赝势对磁系统的描述更准确。然而，对于非自旋极化体系中准芯区电子，NLCC也具有同样的作用。DFT总能量准确表达需要NLCC，如下： Etot?T????Eion????Eee????Exc??? 在赝势的计算式中，电子密度分别来自于芯区电子和价电子。将芯区能量假设为一常数并切不计入计算。用一个价电子密度和由赝势计算得到的离子定域势Eion来代替总电子密度，这样芯区电子与价电子之间所有的相互作用全部转移到赝势上。由此可以推断电子密度线性化只是对动能和简单非线性交换-相关能的一个近似，很明显当芯区电子和价电子在空间很好分离时是一个良好的近似。但如果两个区域电子密度的叠加密切时，计算体系本身就会产生错误，进一步减弱赝势实用性。解决NLCC问题的方法就是调节赝势生成方法以及在固体中计算方法。在产生赝势时每个角动量通道对应一个屏蔽势，并且满足一定的条件，比如规范-保守，赝波函数本征值与全电子波函数本征值相同等。这些屏蔽势（screened potentials）对应的原子赝波函数（atomic pseudowavefunctions）仅表示价电子。从这些波函数可以得到价层赝电子密度（pseudo charge density），通过对势的屏蔽得到“光秃”离子势（bare）： lVion(r)?Vl(r)?Vee[??(r)]?Vxc[??(r)] 由于交换-相关势泛函是电子密度的非线性函数，对自旋极化体系采用这种方法产生的离子势与价电子排列有关。Louie等提出了将上面方程替换为如下表达： 8 CASTEP计算原理---------XBAPRS V(r)?V(r)?Vee[??(r)]?Vxc[??(r)??c(r)] 在屏蔽原子势中减去总交换-相关势。此外，在计算交换-相关势时芯区电荷必须加到价电子中去，这个额外原子状态信息传递给CASTEP，在所有计算中芯区电荷认为（deemed）是相同的，这种做法的一个缺点是在利用赝势计算时芯区电荷很难准确的用Fourier网格表示。而且通常芯区电子密度比价电子密度大，这很容易将与价电子密度有关的影响掩盖掉。以下部分将对部分芯区校正方程建立做介绍，该方法充分的认识到价电子与芯电子密度重叠的区域才是我们感兴趣的。靠近原子核的芯电子密度不会产生物理结果，虽然有如上所述的一些问题。部分NLCC采用一个在特定半径以外与?c一致的函数替代全芯电子密度，在原子核周围这个函数起伏是平滑的。在CASTEP中对一些特定元素在赝势中采用的部分芯区校正使用了数值化的芯区电子密度。在规范保守赝势中虽然有相关的内容，但在计算中并没有采用这个方法。 lionlA Introduction to DFT 第一性原理（The first principle）计算也称为从头算起（ab-initial calculation），由于固体的许多基本的物理性质是由其微观的电子结构决定的，因此通过求解多粒子系统的Schodinger方程，来获取固体全部的微观信息从而预测宏观的性质。利用这个思想建立的能量的哈密顿量非相对论形式可以表示如下： ?2?212?2m??i2?(r,r,r,.......,rrr.........,r)?2???(r,r,r,.......,rrr.........,r)?ke1e2e3ei,n1,n2,njei,n1,n2,njeijmnjjke1e2e3(?????)?k(re1,re2,re3,.......,rei,rn1,rn2,.........,rnj)i1,i2rei?reji,jrei?rnjj1,j2rn?rnj1?j2j1j2i1?i2?Ek?k(re1,re2,re3,.......,rei,rn1,rn2,.........,rnj)考虑到原子核与核外电子质量差别以及电子驰豫时间比原子核驰豫时间要小三个数量级，因此利用 Born-Oppenheimer近似将原子核运动和电子的运动分离，从而将体系波函数划分为电子波函数和原子核波函数两个部分，分别用?和?表示： e2zje2zj1zj2e2?(r,r,r,.......,rrr.........,r)??(r,r,r,.......,r)?(rr.........,r)ke1e2e3ei,n1,n2,nje1e2e3ein1,n2,nj能量的哈密顿量可以分解为如下的两个方程： 9 CASTEP计算原理---------XBAPRS 22e?2zeej2(?2m??i????)?(r,r,r,.......,r?E?(r,r,r,.......,re1e2e3ei)kee1e2e3ei)eii1,i2rei?reji,jrei?rnj i1?i22n?2zje2zj1zj2e21(?2??????)?(rr.........,r)?E?(rr.........,r)jn1,n2,njknn1,n2,njmjnji,jrei?rnjj1,j2r?rnnj1?j2j1j2第一性原理严格求解仅在氢分子中实现了，对于多粒子体系的计算几乎是不可能的。目前均采用不同的 近似方法来实现计算，主要方法有Hartree-Fock近似和DFT近似。在Hartree-Fock近似中体系的哈密顿量表示如下： ETotal?HFi??)???HF???(Jij?Kijiij?iJ ??为第i个电子的Hartree-Fock的轨道能，ij是库仑积分，表示电子静电互斥能，Kij为交换积分。 交换积分所代表的交换能指电子由于自旋平行而引起的电子轨道库仑能量减少的部分。 密度泛函理论(Density Functional Theory)建立了将多电子体系化为单电子方程的理论基础，并且给出了有效势计算方法，是目前研究多粒子体系性质的一种普遍使用的重要方法。 该理论认为对于处于外势场V（r）中相互作用的多电子系统，电子密度分布函数ρ(r)是决定该系统基态物理性质的基本变量。密度泛函理论中体系的能量泛函表示如下： E(?)?T(?)?U(?)?E(?)txc T(?)：Kinetic energy; U(?):classical electrostatic energy; Exc(?):exchange and correlation energy 由上表达式可见体系能量是电子密度的泛函，因此可以进一步将上式表达为： 2?(r)?(r?)?eE[?(r)]??V(r)?(r)dr?T[?(r)]???drdr?E[?(r)]txc2r?r? 在上式中第一项为电子在外场中的势能，第二项为电子的动能，第三项为电子相互之间的库仑能，第四项为交换相关能，最后一项形式是未知的。 系统的电子密度分布可以表达如下： 2??(r)??i(r)i 10 CASTEP计算原理---------XBAPRS 利用上式可以将动能项表示为： T[?(r)]?n??2?i??i2i U(?)表达为： nNN1?z1U[?(r)]????r(r)?i(r)?2???i(r1)?j(r2)?i(r1)?j(r2)???Ra?rr1?r2iaija??aNNza11????(r1)??(r1)?(r2)???Ra?r12r1?r2a??aaVe?r????(r)VN??(r)?VNN2E(?)zaz?Ra?R?zaz?Ra?R? xc形式确定有两种方法：局域密度近似（LDA,Local Density Approximation）和广义梯度近似（GGA, General Gradient Approximation）。在局域密度近似（LDA）中采用了均匀电子气的分布函数推倒出了非均匀电子气的交换-相关能泛函，从而得到Exc(?)的具体形式。从近期计算结果相关报道来看采用局域密度近似（LDA）计算在绝缘体中会产生较大的误差，而且对带隙宽的半导体等得到不正确的结果。采 用局域密度近似（LDA）主要的缺陷现归纳如下： 1． 对光学跃迁带隙预测很差（一般是过低估计带隙宽度）。这虽然对基态性质如电荷密度，总能量以及 力影响不大，但在导带状态计算中却是个大问题，如关于光学性质，运输性质等的计算。在诸如光伏装置等领域的研究中，带隙就是个很重要的问题。采用“剪刀”(Scissors)工具在固体带隙计算中很有用，但对我们未获得实验结果的物质，是不能采用这个方法的。 2． 对类似于二氧化硅这样的电子气分布极不均匀体系，基本假设中关于电子密度分布在空间是缓慢变化 的条件是不满足的，这样的体系采用LDA处理就存在难题。 3． LDA简单的认为计算体系是顺磁性(Paramagnetic)的，对于包含未配对(Unpaired)自旋体系采用局域自 旋密度近似（LSDA）（对自旋向上(spin up)和向下(spin down)的电子分别采用密度泛函计算）是很有用的，比如费米能级(Fermi level)处半填充的系统。 4． 最后一个很少关注的领域就是玻璃陶瓷工业，LDA对弱的结合键（如偶极涨落）很难描述，氢键 (Hydrogen bond)在LDA中也无法获得准确的计算结果。 GGA近似则改进了L(S)DA,将相关交换能确定为电子密度极其梯度的函数，在GGA学派中以Perdew等人认为交换相关能的泛函形式应该以一定的物理规律为基础，构造了著名的PBE泛函。 将电子密度分布函数带入体系能量电子密度泛函中，对泛函变分求极小值，可以得到Kohn-Sham方程： 11 CASTEP计算原理---------XBAPRS ??F?(r)?E?(r)iii?F??? ? 2?V[?(r)]KS ???Exc?(r)V[?(r)]?V(r)??dri??KS????(r)ir?r 交换-相关能可以按照下式计算： ?E?xc[?]???(r)EE[?(r)]xc[?(r)]dr :number of particles; xc:exchange-correlation energy per particles in an uniform electron gas ; ?(r)?Exc :distribution function of electron density. ???(r)称为交换-相关势和，表示为： ?xc[?(r)]????(r)?Exc 在Castep计算中采用了周期性边界条件，单电子的轨道波函数满足Bloch定理，采用平面波展开式有： ?i(r)?eiK?R?i(r) 周期性边界条件下的波函数扩展为一系列分离的平面波波矢，这些波矢与晶体的倒易点阵矢量相联系。 ?i(r)??Ci,Ge2.2 晶体光学性质的计算基于以下原理： 12 iG?R CASTEP计算原理---------XBAPRS 电磁波在真空以及某种材料介质中传播时差别可以用一个复数式的折射指数来表示： N?n?ik2k???c..............221?NR?1?N 在真空中N为实数，而且其大小为1；在其他介质中时若材料对于光是透明的则是一个纯实数，虚部对应材料的吸收系数（Adsorption Coefficient）。它们之间的关系方程2所示： 吸收系数表示的是电磁波通过单位厚度的材料时能量的衰减分数，通常可以用材料焦耳热的产生来衡量。 反射系数（Reflection Coefficient）可以简单通过将垂直光束照射材料的表面引起 (n?1)2?k2?......................322(n?1)?k 在计算光学性质时一般先计算虚部的介电常数，其他的性质与介电常数之间建立关系。虚部介电常数计 算式由下方程确定： ??n?k??N......................512这样折射指数的实部和虚部以及介电常数之间的关系可以写为： 222 ?1?n2?k?22?2nk...............6 光导率（Optical conductivity）也是一个普遍用来描述材料光学性质的物理量。光导率的表达式为方程7： ???1?i?.....................72 这个参数用来描述金属的光学性质，但在CASTEP中将计算范围扩大到了绝缘体和半导体。计算过程的主要的区别在于前者的光学谱中IR部分与内部能带之间的转变密切相关，而者则在计算内时并没有完全考虑到这些因素。 从虚部介电常数可以进一步得到材料电子的能量损失函数(Energy Loss Function)，它描述了电子通过均匀的电介质时能量的损失情况，计算式如下所示： ?1Im(?(?)).....................8 在实验中我们可以测定的光学性质参数有吸收系数?(?)和反射系数R(?)。从理论上而言，得到这些 13 CASTEP计算原理---------XBAPRS 参数以后可以将方程2、3、4表示为复数的形式之后得到表达式1中的实数部和虚数部。但在实际情况下由于入射光源的复杂性，而且晶体结构中极化效应使得材料介电常数并非是各向同性的。此外材料表面几何结构也不是理想的平滑表面。这些因素就限制了对其光学参数的预测。在CASTEP中提供的光学性质的计算支持体系极化，但状态只能在同种自旋间相互转化。 晶体中声子和电子之间的相互作用可以用电子基态波函数中包含的含时微扰项来表示，声子电场扰动引起了电子函数占据态和未占据态间的转变（磁场引起的效应要弱一个因数V/C），这些激发态（激子）聚集态称为等离波子。单独的态激发称为单粒子激子，这些激子对光谱产生的结果是导带和价带的状态密度之间的连接可以通过选择合适的加权性矩阵元素来实现。在CASTEP虚部介电常数的计算按照方程9进行： 22?2e?cu,r???(Ec?E??E).....................9??(q?O?,??)????kkkk20k,?,cu u矢量定义光束电场的极化性质。这个表达类似于含时微扰的Fermi-Golden定理，?2(?)可看作真实 占据态与未占据态之间转换的细节。介电常数就描述了一种因果效应，它的实数部和虚数部之间由Kramers-Kronig变换相联系。利用这个变换就可以得到介电常数的实数部?1(?)。 用于描述电子态转变的位置算符矩阵元素通常用动量算符矩阵元素来表示，这样可以在倒易点阵空间直接的进行计算。局域势函数会影响计算，在CASTEP计算中一般采用非定域势函数。本文在进行BFGS晶体结构几何优化时就选择了非局域势函数。经过矫正后的矩阵元素可以描述如下： 11??cP???c(V,r)??...........................10r?c???i?m??kkkkknlk 利用超软赝势（Ultra soft Pseudopotential）计算时会增加额外矩阵元素，在目前CASTEP计算中这部 分矩阵元素并没有涉及。采用规范保守势计算结果发现与采用超软赝势计算符合的很好，因此额外的那部分矩阵元素对于计算结果的影响不大。 晶体光学性质IR部分受能带内部的影响较大，采用经验Drude表达形式就可以精确地描述这个影响。 ?0?D(?)?i?1??D Drude校正的光导率?0和Drude限制系数?D与材料许多实际参数有关，一般这些参数可以通过实验得到。结合上式和式7就可以了解介电函数中Drude的贡献，同样可以得到在其它光学常数中的分布。Drude限制参数描述了计算过程中未涉及因素引起光谱宽化现象，比如电子间的散射效应（包括Auger效应）、电子与声子之间的散射效应以及电子与晶体结构缺陷之间的散射效应等。 在CASTEP中光学性质计算结果的准确性与下列因素有关： 1.导带数量（Number of conduction bands）：直接决定了Kramers-Kronig变换的准确性。 2.截止能量(Energy cutoff)：体系能量进行迭代计算过程中，电子基态能量本征值精度直接影响能带结构以及光学性质，提高截止能量的数值可以提高计算精度，可以得到更准确未占据态的自恰电荷密度和震 14 CASTEP计算原理---------XBAPRS 动自由度。 3.迭代计算中K点数量(Number of k-points in the SCF calculation)：与截止能量对体系基态能量计算影响一样，K点数量越多，迭代计算能量越准确。 4.积分Brillouin zone K点数量(Number of k-points for Brillouin zone integration)：在计算光学性质矩阵元素时Brillouin zone选取的K点数量应当是合适的，与电子能量相比，矩阵元素在Brillouin zone变化更快，因此必须选取足够数量的K点来提高矩阵元素计算结果的准确性。 从目前计算结果对比来看，提高上述参数的准确性时，光谱中特征峰可以快速地达到实际的要求。当然CASTEP中对光学性质的计算还有不少的局限性，电介质极化引起的局域场效应在现在计算中被忽略了，这对光谱计算有一定的影响，但在目前计算方式下将是无法进行的。准粒子和DFT能带带隙以及激子等都会影响计算结果。 状态密度在Brillouin zone区的表示： 给定能带n对应的状态密度Nn(E)定义为： dkNn(k)??3?(E?En(k)) 4?En(k)描述了特定的能带分布情况，积分在整个Brillouin zone进行。另外一种表示状态密度的方法基于 Nn (E)dE与第N级能带在能量E到E+dE范围内允许波矢量数成比例。总体状态密度N（E）就是对所有的能带允许电子波矢量求和，从能带极小值积分到费米能级就得到了晶体中包含的所有的电子数。在自旋极化体系中状态密度可以用向上自旋（多数自旋(majority spin)）和向下自旋（少数自旋(minority spin)）分别进行计算，他们的和就是整体状态密度分布，它们的差值称为自旋状态密度分布。借助于状态密度这个数学概念可以直接对电子能量分布进行积分而避免了对整个Brillouin zone积分。状态密度分布经常用于快速直观的分析晶体的电子能带结构，比如价带宽度、绝缘体中能隙以及主要特征谱峰强度分析，这对于解释实验各种谱数据有很大的帮助。状态密度还可以了解当晶体外部环境如压力等发生变化时电子能带的变化情况。 状态密度数值化计算方法很多，最简单的方法是对各个能带电子能级进行采用柱状图取样Gaussian拟和。用这种方法绘制的状态密度分布图不存在类似于van-Hove奇点尖锐分布，但只需要少量的K点即可。其他的准确方法基于对Brillouin zone参考点之间采用线形或二次方内叉法。目前最可靠和普遍使用的方法是四面体叉入法，但这种方法与Brillouin zone网格特殊点是不融合的。因此CASTEP使用了由Ackland发展的简单的线性内叉法，对Monkhorst-Pack倒易基组平行六面体采用线性内叉法，能带能量组合基组进行柱状取样。 2.4 偏态密度（PDOS）和局域状态密度(LDOS) 偏态密度（PDOS）和局域状态密度是一种分析电子能带结构有效的半经验方法。局域状态密度表示了体系中不同原子在各个能谱范围内电子状态分布情况。偏态密度（PDOS）进一步将上述分布以角动量贡献进行量化分析。了解状态密度分布峰值中S、P和D轨道贡献是很有用的。LDOS和PDOS提供了一种定量分析电子杂化状态的方法，对于解释XPS和光谱峰值的起源很有帮助。PDOS计算基于Mulliken population分析，每个给定原子轨道在能带各个能量范围内分布均表示出来，特定原子所有轨道的状态密度分布和以LDOS表示出来。与整体态密度计算相似，采用了高斯混合算法或线形内叉法。 Brillouin zone积分取样 大快固体中电子状态只允许存在于由边界条件确定一系列ｋ点中，固体周期性结构中包含了无限数量的电子，这对应于无限数量的ｋ点。无限数目的电子波函数计算利用Bloch定理转变为用有限数量ｋ点计算有限数量的波函数。每个ｋ点处电子占据态都会对电子势有贡献，因此在理论上要进行无限数量的计算。对于十分临近的ｋ点，它们的电子波函数几乎是完全相同的，因此在ＤＦＴ表达中对所有ｋ点求 15 CASTEP计算原理---------XBAPRS 和（等价于对整个Brillouin zone积分）可以采用有效的离散化数值计算，即在Brillouin zone选取有限数量的特殊点。进一步考虑到对称性，只对Brillouin zone无法简并的部分才计入计算过程。Payne以及Srivastava and Weaire等人的文献提供特殊ｋ点选择方法以及求和加权的评论。采用上述方法以后，选用很少的ｋ点对绝缘体电子状态计算就可以获得对电子势和总能量准确的近似。对于金属体系而言为了得到费米能级准确性，需要更致密的ｋ点数量。采用更多ｋ点数量就可以减小因Ｋ点数量限制而产生的对总能量计算的误差，与获得基组数量方程收敛方法类似。当对对称性不同的两个体系的能量进行对比时，与ｋ点取样相关的计算收敛精度要更高，例如比较ＦＣＣ或ＨＣＰ结构相对稳定性。在这种情况下计算误差是不可避免的，因此能量必须达到绝对收敛精度。 要注意的是，体系总能量不会因ｋ点数量的不同而发生变化，因此即使收敛精度很低时能量计算也一样，这就与平面波基组截止能量的收敛计算不同，后者平面基组增大时总能量会减少。 Monkhorst-Pack特殊点（ special points） Monkhorst-Pack发展了一种目前普遍采用的特殊ｋ点产生方法，最初只在立方体系中使用，后来 Monkhorst-Pack将其进一步扩展到了六方晶格中，在倒易空间沿着坐标轴生成均匀规则分布的ｋ点网络。Monkhorst-Pack网络采用三个积分来定义，qi where i=1,2,3，确定了与主坐标轴之间的偏差。这些积分得到了下面的一些数字： ur=(2r-qi-1)/2qi where r varies from 1 to qi. The Monkhorst-Pack grid is obtained from these sequences by: kprs=upb1 + urb2 + usb3 q1q2q3这个基组不同点进一步调和，对调和基组中的特定点按照其镜像对称点进行加权性取样。 在对基组中所有点调和前，可以增加一个常数变化，应用于六方点阵结构时，在沿a and b 轴方向所有点产生一个轻微修正的结果。 up=(p-1)/qi Where p varies from 1 to qi. 计算材料学报告中应当注意的问题： 随着新一带材料学计算软件的不断开发和更新，采用计算机来模拟和预测材料的性能已经成为计算材料科学中的前沿热点，每年全世界有数百篇与此相关的论文发表。但这些模拟的结果很大一部分无法得到很好的再现，因而存在大量的自相矛盾的信息。在这里实际上很难判断在某一次计算中采用的模型，算法是否是存在问题的，Ann E Mattsson1, Peter A Schultz等人提出了如何才能获得有意义的模拟结果，从计算方法，平面波基组，能量截止，赝势函数，与计算性质相关的超晶胞结构的建立以及周期性边界条件的设定等一系列的问题都对最终的计算结果产生影响，因此当论文中出现的结果出现矛盾时就需要通过对计算细节的描述来判断其正确性。一般而言计算结果是冗长的，因此有必要将其与相应的论文在网络上发表，利用因特网来让研究人员能够获得这些细节信息，从而对论文的计算结果进行重复和验证。为此，他们提出以下的指导性意见： 影响计算结果精度的因素： 1.赝势选用（ PPs）: If used, identify them. Any deviation from standard, published PPs should be described in 16 CASTEP计算原理---------XBAPRS sufficient detail for the work to be reproducible. 2. k points: Report the sampling that was used and which convergence tests were performed. 3. Basis sets/kinetic energy cutoff: Basis sets should be identified and, if non-standard, the reason for using them given. If feasible, a calculation should be repeated with another appropriate basis set and the result reported. If plane waves are used the cutoff should be given along with the convergence it affords. 4. Trajectory length and time step in AIMD: Motivate why they are appropriate. 5. Equilibration in AIMD: How are the initial configurations prepared? 6. Fictitious electron mass in CPMD: Report and motivate the choosen mass. (b) Factors affecting accuracy: 1. Functional: First and foremost, which functional was used? It is a good practice to repeat some of the key calculations using a second functional and report the result. This provides an estimate of the accuracy as well as information highly important Topical Review R27 for further development of functionals. Even negative results are valuable—which functional not to use for a specific system. 2. System size: Is the property under study converged with respect to the size of the super cell or cluster? 3. Relaxation: Report on whether only volume or full relaxation was used. Justify the reliability of any results obtained for an unrelaxed system. 4. Boundary conditions: The exact treatment of the boundary conditions should be described for a simulation where the system is not a simple bulk crystal. 利用CASTEP模拟计算实例 一，计算本征半导体硅的能带结构和状态密度等性质 计算过程分为三个步骤：首先是建立硅的晶体结构计算模型，这个可以在MS物质结构数据库中调用即可。在计算时为了节省时间，减少计算量将硅的普通的晶体转化为原胞结构，一个原胞中包含9个原子。节下来是对晶体原胞结构进行几何结构优化，当然其中也含盖了对体系总能量的最小化。结构优化过程中的两个图表文档分别表示了优化步骤中体系能量的变化和收敛精度，判断收敛是否成功就要查看最终完成计算后，能量的收敛精度是否达到了事前的设定值。最后是计算性质，在计算状态密度时可以计算不同原子各个轨道按照角动量分布的偏态密度（PDOS），当体系是自旋极化时，偏态密度（PDOS）中包含了体系多数自旋(majority spin)和少数自旋(minority spin)的偏态密度（PDOS）。光学性质的计算是模拟中的一个难点，从目前发表的文献来看，影响光学性质计算的因素很多（见光学计算原理部分，对此有详细描述），在研究体系有充足实验数据的条件下，可以对能带采用“剪刀”的工具对能带带隙进行刚性的调整，获得与实验结果符合较好的结论。但对于初学者而言，这个工具一般是不推荐使用的。作者对于硅的计算完全按照上述方案完成。详细的计算结果和计算方法见本文所附带的专门文章。 二，搀杂半导体InP性质计算 第三主族和第五主族元素之间形成的半导体，目前越来越受到的重视，在纳米材料中，各种纳米电子器件如场效应晶体管，半导体纳米量子阱，纳米量子点激光器等均广泛采用了诸如AlAS InP等材料，本文对InP能带结构、状态密度以及光学性质进行了计算。计算步骤与前文描述相同。详细结果见文章二。 三，FeS2性质计算 二硫化亚铁是一种受到广泛研究的窄带隙的半导体，其能带带隙为0.95eV。肖奇等人也采用CASTEP对二硫化亚铁整体状态密度和（100）晶面双层超结构状态密度的计算结果进行了对比，发现了表面态对状态密度峰的分裂。作者也首先建立了二硫化亚铁的晶体结构，对优化后的结构也进行了计算，得到能带带隙的较准确结果，但在能带的顶层出现了文献中未出现的新结构，因此还需要其他文献进行证实。作者也建立了双层的二硫化亚铁（100）晶面的超晶胞结构，但限于计算能力，只对结构进行了分子力场的初步结构优化。详细结果见文章三。 四，三氧化二铝性质计算 17 CASTEP计算原理---------XBAPRS 三氧化二铝是广泛用于复合材料中的一种附加材料，在电子工业中用于衬底材料。作为陶瓷原料更是普遍。三氧化二铝有多种晶体类型，目前广泛得到研究的是α-三氧化二铝，作者计算了它的能带结构和状态密度分布以及电子密度分布情况，计算结果与实验结果相比较是可靠的，电子密度分布揭示了化学键的性质。详细结果见文章四。 五，其他几种半导体材料能带结构的计算 作者也计算了几种目前普遍使用的半导体材料的能带结构，晶体结构在计算前是经过了结构优化的，某些计算的能带带隙并不理想，与实验数值相比较，差距较大。但发现了能带结构和计算晶体结构特别是化学键类型间的关系是密切的。 通过对于几种类型的半导体能带结构和状态密度的计算表明他们与原子轨道杂化类型，原子间成键类型等均有关系，计算几种半导体分别是：本征半导体Si,离子型窄带隙半导体ZnO and Cu2O，搀杂n型半导体BN。 从杂化轨道类型来看，硅为sp3杂化，BN为sp2杂化。其余两种是离子型晶体，化学键主要成分是离子键。从能带结构分析，离子型半导体费米能级以下的部分能带形状平滑，而共价键杂化类型的半导体在费米能级以下部分为抛物线型。在费米能级以上的部分两者差别不大，均为抛物线型。状态密度(DOS)图来看，ZnO and Cu2O型半导体各个分波区分是很明显的， 杂化型半导体状态密度各个分波区分并不明显，一般为连续型，原子轨道混合在这些半导体中是很明显的。 18 CASTEP计算原理---------XBAPRS 六，MnN和MnAs自旋状态密度分布与晶体结构常数间的关系 R. de Paiva1, J. L. A. Alves等人文献中，研究了闪锌矿结构的MnN 和MnAs自旋密度分布随晶体结构常数变化的关系，上述两种物质的闪锌矿结构并不是它们的稳定结构，但在这种结构中Mn以四面体配位性质在许多二元磁性材料均有体现，人们相信，正是锰元素这种配位环境形成了独特的磁性质。文献作者采用第一性原理的DFT理论分别用GGA和LDA分别计算了随着晶体尺寸变化，自旋密度分布的变化情况，他们发现上面两种物质自旋密度状态分布结构与晶体尺寸密切相关，当MnN尺寸大于0.490nm,MnAs大于0.571nm 延伸闪锌矿结构是半金属型的,即多数自旋(majority spin)密度分布在费米面处是连续的,少数自旋(minority spin)密度则在费米面处是绝缘体型的。他们计算中晶体平衡尺寸分别是： MnN：ao = 4.19 A° (LDA), and 4.30 A° (GGA)); MnAs： the magnetic moments are 2:51B and 4:01B respectively at the equilibrium lattice parameters ao = 5:32 ° A(LDA) and 5.71 °A (GGA). 作者计算均采用GGA，MnN：0.425nm；MnAs：0.5977nm。自旋密度分布和能带结构基本与文献结果一致，在MnN中，自旋密度在费米能级附近随着尺寸变化是很明显的，但MnAs中则与文献结果不一致。 19 CASTEP计算原理---------XBAPRS 现将文献结果与作者结果列图如下： 从多数自旋(majority spin)密度和少数自旋(minority spin)密度图中可以明显的得到体系电子电导的单自旋极化现象，半金属特性是 明显的。但MnAs计算结果与文献差别较大，特别是自旋密度与晶体随尺寸变化的关系不明显，作者计算 20 CASTEP计算原理---------XBAPRS 了在0.571nm,0.569nm and 0.550nm的情况下的状态密度情况，仍然没有发现自旋密度分布在费米面处明显的变化。特别是文献中提到的少数自旋(minority spin)密度在晶体尺寸小于0.571nm时会在费米面处产生能隙。 七，孪晶界（twin boundary)的建立 作者仿照CASTEP文献提供的帮助资料建立了硅的一个（310）的孪晶界（twin boundary)，主要步骤有：晶面切割，选择合适的切片厚度，按照建立表面模型的方法建立（310）和（-310）孪晶界（twin boundary)结构，并且采用分子力场进行优化。如下所示： 21 CASTEP计算原理---------XBAPRS 八，DPD动力学模拟双相分离 作者利用这一模块初步进行了双相分离动力学模拟，利用这个模块模拟了水和油两相的分离情况。 关于LDA和GGA的比较 这里主要讨论LDA和GGA在计算不同的模型时的结果比较，可以对以后计算中选用较合适的泛函形式提供一定的参考。下面文字由作者摘录并且翻译自Ann E Mattsson1, Peter A Schultz等人的文章。 John Perdew建立了目前广泛使用的泛函，他将Jacob’s ladder作为泛函发展过程中的坐标，用这个概念来对泛函分类是很有用的。每上一个阶梯都会增加泛函形式的复杂性。 LDA局域密度近似（LDA）：局域密度近似（LDA）是第一阶梯。它仅仅采用空间点r处的电子密度n（r）来决定那点交换-相关能密度的形式。交换-相关能密度由密度相同的均匀电子气完全确定。泛函的交换部 22 CASTEP计算原理---------XBAPRS 分就准确的用均匀电子气的微分表达。各种不同的局域密度近似（LDA）仅仅是相关部分表示方法不同，所有现代应用的局域密度泛函都基于Ceperly和Alder`s在８０年代对均匀电子气总能量的Monte Carlo模拟。Perdew-Zunger(PZ),Perdew-Wang(PW) and Vosko-Wilk-Nusair(VWN)对ＣＡ数据进行了不同程度的拟和。 广义梯度近似（ＧＧＡ）：ＧＧＡ是Jacob阶梯的第二个台阶，将电子密度的梯度也作为一个独立的变量（|?n(r)|），在描述交换－相关能方面，梯度引入了非定域性。ＧＧＡ泛函包含了两个主要的方向：一个称为“无参数”，泛函中新的参数通过已知形式中参数或在其它准确理论帮助下得到。另外一个就是经验方法，未知参数来自于对实验数据的拟和或通过对原子和分子性质准确的计算。Perdew,Burke and Emzerhof(PBE)以及Perdew-Wang from 1991(PW91)是无参数的，在量子化学中广泛采用的GGA，比如Becke,Lee,Parr and Yang(BLYP)是经验性。LYP校正采用了密度的二阶Laplace算符，因此严格上讲属于Jacob阶梯的第三阶，但通常仍然归类为GGA。 Meta-GGA:第三阶梯的泛函使用了电子密度拉普拉斯算符或动能密度作为额外的自由度。Tao,Perdew 及其同事最近构建了无参数的meta-GGA，TPSS泛函。经验的meta-GGA采用了通过拟和得到的相应参数。 Hyper-GGA：增加精确交换(exact exchange)构成了Jacob台阶的第四阶梯。每个粒子的能量可以通过SE多体波函数计算得到，这就是Hartree-Fock(HF)交换公式。在上述公式描述中采用KS单粒子本征函数而不是HF轨道，通常就叫精确交换。 广义随机相近似(Generalized random phase approximation):泛函的最高阶梯采用了EXX和精确部分交换。 采用其他的泛函分类方法不如Jacob阶梯明晰。一种方法是亚体系泛函，体系不同的部分采用了经过不同修饰的泛函来描述每个部分主要的物理相关作用。另一种方法是将分子间相互作用(Van der Waals)描述也含盖到泛函中。表面校正方法的发展和引入可看作最原始的亚体系泛函，精确的表面能可以从“jellium”表面获得（一系列表面模型），LDA，GGA或这些体系中的其他的泛函表面能的误差都可以计算。通过假定一个特定的泛函可以给出在真实表面和与之相关的“jellium”模型相同的表面能，对真实体系DFT计算的表面能采用“尾部”校正来进行评估。 通过对实验数据拟和或SE计算的准确值建立的泛函是另一种称为杂化(hybrids)泛函。在杂化密度泛函中，一个混合的HF交换和GGA交换使用时与GGA或（和）LDA交换校正相结合。通过对大量实验分子数据的拟和来决定混合程度。在量子化学中使用最成功的泛函，B3LYP（Becke三参数杂化泛函与LYP泛函相结合方法）就是杂化形式。杂化TPSSh是TPSS与HF交换的混合形式，杂化泛函和其他泛函使用性能的不确定性与拟和时采用的原子和分子种类以及性质有关。各种泛函近似的准确性会随着研究材料体系和性质而发生变化，对所有研究体系而言，没有一个泛函是最好的，甚至谈不上“好”，在物理与化学领域所观察到的结果反映了各种泛函基本的不同之处。与物理和化学有关的领域选用合适的泛函是很有挑战性的问题，比如表面的化学反应。以下部分将讨论与泛函近似有关的焦点问题。 选择合适泛函计算：某个计算模拟中泛函的选择过程是一个应当引起足够重视的环节，在同一个软件包中同时的使用各种可能的泛函和仅发表与已知实验结论一致的计算结果违背了“预测”计算前提，实际上这种一致性在很多时候仅仅是一个巧合。 在泛函选择上没有足够可以保证的条件，但在下面我们将给出一些指导： 众所周知的是爬上泛函的Jacob阶梯可以系统的提高计算结果，对多分子体系采用GGA在多数情况下胜过LDA。目前已经广泛的认识到对水的描述采用PBE和BLYP是优先选择，而不是LDA。然而，不总是这样。对某些体系和性质的计算采用LDA要比GGA好，特别是表面能和许多氧化物性质的计算。因此，在计算中采用不同形式的泛函同时计算 来获得较为准确的评估。在任何的DFT计算中所采用的泛函形式都应该以书面形式报告，在实际应用中一个LDA的准确度很少受与所选用的相关泛函影响，在LDA计算中采用的校正泛函必须说明。在1980年以前采用的LDA相关泛函并不是基于对CA电子气计算建立的，因此与现在采用的泛函相比计算结果的差别是比较大的。PBE和PW91泛函采用了类似的方法来构建，计算结果的差别也不大。与PW91相比PBE在芯区有效势更为平滑，后者数值上不稳定。一般对这些并不进行区分， 23 CASTEP计算原理---------XBAPRS 通记作GGA，不过我们和其他的研究人员鼓励还是应当对计算时采用的PBE 或PW91予以区分，不为人知的事实是一些GGA相关泛函，比如PBE和PW91，是在PZ，PW和VWN中嵌入一个LDA相关泛函，如果可能的话也应当予以说明。密度泛函理论已经广泛的在各种计算中得到应用，已经有很多的文献阐述了在各种不同类别的材料研究中采用不同泛函计算的问题。就像其他研究一样，论文提供了对某个新问题成功的选用泛函的详细情况。但这些并不都十分可靠，比如说铝空位形成能的计算。GGA（PW91）计算可以得到关于体系性质较好的结果，如点阵常数、弹性模量，结合能等，因此GGA对铝的计算中优于LDA。但对单个空位，这种假设是不成立的。LDA计算的单个空位的形成能，0.7eV,与实验测定值0.67(3)eV相比要好于PW91计算值0.54eV。事实是单个空位形成能接近于表面性质而不是体系性质，而如上述LDA在计算表面能方面比GGA要好。当然对单个空位形成能进行一个表面校正，GGA性能就会超过LDA。 铝的这个例子就说明了在任何计算中研究使用不同泛函的价值所在。而且，在计算中如果选用的泛函并没有给出满意的计算结果，我们并不能因此轻易的放弃它。将这些否定的结果与计算成功的泛函一同发表对以后DFT研究相同或类似计算提供了有价值的信息。对同一个计算过程采用不同的泛函可以提供计算结果理论误差的界限。对于某些问题，特别是只对体系量化分析时，采用不同的泛函可以给出可对比的结果，误差范围可以小到不会改变体系的性质。对在其他方面的应用，特别是要得到量化信息，误差界限会增大使得独立的预测就无法给出。然而，即使利用不同的泛函进行的DFT的计算也没有获得结论，额外的信息特别是实验数据就可以确定那个泛函是可靠的并且得到目前研究体系的更好的信息。合适的泛函在材料计算模拟中准确性是可信的，但并不总是这样。现在用于计算的泛函对表面性质的计算并不可靠，对分子间相互作用或色散作用这一点是很重要的。另一个问题就是LDA或GGA泛函对半导体或绝缘体能带带隙计算结果的可重复性不强。带隙问题从大范围上来说是“自相互作用”误差引起的，采用EXX来有效的消除“自相互作用”误差极大的激励了这个领域的发展。然而就EXX本身还不够，引入一个融洽的校正泛函就完成全相关-交换泛函的建立。到目前为止我们还没有涉及到DFT计算中自旋问题，多次经验表明对于绝缘体中的缺陷选用自旋极化势是得到满意结果的关键条件。作为一种主要的改进方法，仅仅采用自旋极化也是不够的。两个最近的研究事例：一个是石英中铝的置换，另一个是氯化钠表面的缺陷研究，都说明现在的泛函自旋极化密度过于非定域化。“自相互作用”误差又一次成为了一个计算缺陷，不过采用EXX可以得到改善。我们的最后一个例子是物理和化学交叉的领域---在大块材料表面分子的吸附研究。实验中发现一氧化碳分子会直接吸附在单个的Pt（111）表面上。当前泛函LDA或GGA以及计算方法，全电子或赝势都预测一氧化碳分子会驻留在一个空位处与三个Pt原子配位。各种各样不同的解释都用来描述理论预测与实验之间的差别，对目前泛函在这个问题上的失败有待于未来获得解决。这就说明利用现在的泛函来进行DFT计算，即使在最理想的计划和控制下，也并不能总是获得正确的结果 最近作者计算了WC性质，获得了与文献完全一致的结果。 几种常见的赝势形式及表达式 24 CASTEP计算原理---------XBAPRS By XBAPRS All Rights Received. For more informations ,please correspond to xbaprs 2005-8-24 25 CASTEP计算原理---------XBAPRS 分子力场计算方法概论 分子力场计算等级介于第一原理计算和经典分子动力学之间,是一种半经验的算法,在MS中Forcesite就是这种算法的模块。与经典分子动力学相比，半经验算法准确性似乎更好一些，因为半经验算法中原子之间相互作用能参数是通过严格的第一原理计算得到的，可能采用的是Hartree-Fock的波函数算法，也可能是基于DFT的算法。在经典分子动力学模拟中采用的一般是各项同性的对势函数，比如Lennard-Jones势，还有一些采用指数衰减来模拟化学建的势函数如Morse势等。当然最近在考虑到电子结合相关作用的基础上又引入了对泛函势，二体相互作用也可以扩展到多体相互作用，如团簇势。分子力场的算法实际与经典分子动力学中对原子相互作用的描述有些类似。对体系总能量的计算而言，分子力场中由以下几个部分组成，即键能(Valence)，非键作用能(Non-bond)和交叉相互作用能组成(Cross terms) Etotal?Evalence?Enon?bond?Ecrossterm Cross-terms主要是组成键能各部分之间的耦合作用。Valence这部分能量还可以进一步划分为下面五种能量： Ebond：Bond stretch; Eangle：Valence angle bending; Etorsion：Dihedral angle torsion; Eoop：Inversion (out of plane interactions); EUB：A Urey-Bradely interaction. Ebond?Eangle?Etorsion?Eoop?EUB EUB是一种间接的三体作用，一般是存在于同时连接在同一个共用原子两端的原子之间。Eoop是共价键 特有的相互作用。 Cross-terms是为了提高计算准确性引入的校正项，存在这些项的分子力场计算方法一般是现代力场算法，这些校正项用来表示邻近原子对化学建以及键角的畸变。校正后的分子力场算法可以更好的获得与实验数据相符合的结果，比如分子震动频率谱。交叉相互作用项主要有一下几部分：Stretch-Stretch、Stretch-Bend-Stretch、Bend-Bend、Torsion-Stretch、Torsion-Bend-Bend、Bend-Torsion-Bend and Stretch-Torsion-Stretch。计算经验表明这些耦合作用项的存在对计算结果不一定有利，对于畸变很大的晶体结构或分子结构，所得到的优化结构是不真实的，因此对这类计算还是选择初级的分子力场算法，防止体系陷入局域最小状态。 Non-bond作用有三部分构成：Van der Waals (VdW)、Electrostatic (Coulomb) and Hydrogen bond (h-bond)。 Enon?bond?EVdW?ECoulomb?EH?bond 分子力场计算中一般的表达式如下所示：计算式中前四项分别表示:Stretch bond (b), Bend angles (?) away from their reference values, Rotation torsion angles (?) by twisting atoms about the bond axis that determines the torsion angle, Distort planar atoms out-of plane formed by the atoms they are bended to (?)。 26 CASTEP计算原理---------XBAPRS 后面的五项就是Cross-terms，表示前面四种作用之间的相互耦合作用。最后一项是非键作用，包括Lennard-Jones势形式相互作用以及Coulomb形式的静电作用。 V(R)??Db[1?exp(?a(b?b0))]2??H?(???0)2??H?[1?Scos(n?)]b??`??H??2???Fbb`(b?b0)(b`?b0)???F??`(???0)(?`??0`)?bb`??`???Fb?(b?b0)(???0)??F??`?(???0)(?`??0`)cos????F??`??`b???`??` ???[ij?iAijr12ij?Bijr6ij?qiqjrij]分子力场计算方法特点 分子力场计算方法与严格的量子力学算法相比，最大的特点就是处理体系可以很大。在分子力场的算法中很多的参数来自于广泛实验数据的拟和和严格量子力学计算，因此这类算法不仅在重复实验结构方面很适用，而且部分量子效应也包含在这些半经验的算法中了。分子力场算法特点可归纳为以下几点： 1． 处理体系大小与量子力学计算相比要大几个数量级，因此采用分子力场算法可以模拟凝聚态分子体 系，生物大分子（蛋白质，多肽结构），晶体结构，有机和无机化合物分子等，计算方面主要是计算一些与微观量子效应不敏感的性质，如相组织的分离，状态方程以及键能等。 2． 可以对构成体系总能量的各个部分进行单独分析，比如前面提到的能量计算中键能由成键作用，非键 作用以及交叉作用等构成，可以单独分析这些作用产生的物理效应。 3． 在分子力场计算中可以方便的引入各种约束条件，比如固定结构中部分原子的坐标，固定晶体晶格常 数，晶体体积等，这些约束在研究类似于界面问题时是很重要的。 分子力场中无法计算的体系有： Electron transition state (Phonon adsorption)； Electron transport phenomena； Proton transfer (acid/base reactions). MS中分子力场支持的算法简介 MS中支持分子力场算法的模块主要是Discover 和Forcesite这两个，其中Forcesite提供了目前几乎全部的分子力场算法，比如COMPASS、PCFF、CVFF以及Universal等。现对这些分子立场计算方法适用范围做介绍： 1．Consistent force field-COMPASS and PCFF 自洽场分子力场计算方法包括CFF91,PCFF,CFF以及COMPASS，是第二代半经验力场算法，其能量计算包含了交叉耦合作用校正。这类计算方法中的参数主要由以下物质实验或计算数据得到：含H,C,N,O,S以及P的化合物，主要是由这些元素相互之间形成的物质；卤素原子以及离子；碱金属阳离子以及某些在生物化学领域内重要的金属二价阳离子。PCFF基于CFF91构建，主要增加了聚合物，金属以及分子筛体系方面的计算参数。COMPASS是目前最新一代的分子力场计算方法，其开发初衷就是用于计算凝聚态体系性质的，在金属氧化物以及传统的化合物计算方面已经取得了很大的进步。总的来说自洽场分子力场计算方法还是适合于计算主族元素之间形成的化合物，而且这些物质成键必须是正常的方式。 2．CVFF and CVFF auq CVFF是一种典型的较早应用的分子力场算法，在MS的Discover模块中在进行分子动力学模拟时默认计算方法就是CVFF，这个算法包含了部分非简谐作用以及能量耦合作用项，但不是全部。CVFF算法应用广泛，因此其计算是可靠的，主要用于计算多肽和蛋白质结构。CVFF auq是CVFF的增强版本，增加了非键 27 CASTEP计算原理---------XBAPRS 结合的校正，计算范围也扩大到硅酸盐，铝硅酸盐和粘土等。 3．Dreiding 这是一类比较特殊的分子力场计算方法，适用范围也很广，可用于计算有机大分子，各种有机化合物以及由主族元素构成的无机化合物，同时还有一个特点就是可以计算一些缺乏实验数据的新结构化合物，包括成键特殊的物质，但必须是以上物质范围内的。 4． Universal 适用范围最广泛的计算方法，对元素周期表内的所有元素适用，因此一般结构都可以采用这个算法。 28