能带结构和态密度图的绘制及初步分析
前几天在QQ的中和大家聊天的时候，发现大家对能带结构和态密度比较感兴趣，我做计算已经有一年半了，有一些经验，这里写出来供大家参考参考，希望能够对初学者有所帮助，另外写的这些内容也不可能全都正确，只希望通过表达出来和大家进行交流，共同提高。
MS这个软件的功能确实是比较强，但是也有一些地方不尽如人意的地方。（也可能是我对一些结果不会分析所致，有些暂时不能解决的问题在最后一部分提出，希望大家来研究
研究，看看有没有实现的可能性）。
能带结构、态密度和布居分析是很重要的内容，在
分析能带结构和态密度的时候，往往是先作图，然后分
析。
软件本身提供的作图功能并不是很强，比如说能带结构
（只能带只能做point图和line图），不美观不说，对于
每一个能带的走势也不好观察，感觉无从下手。所以我
一般用origin作图（右图是用origin做的能带图）。能带
结构和态密度的作图过程请参考我给大家提供的动画。
接下来我们先开看看能带结构的分析和制作！
第一部分：能带结构
这个部分打算先简单的介绍一下能带的基础知识，希望能对大家有所帮助，如果对能带了解比较深入的朋友，可以跳过这个部分内容，之中不当之处请勿见笑。^_^
第一个问题是：
1、能带是怎样形成——轨道和一维体系的能带。
这是最基本的一个问题，我们要对能带结构进行分析，首先要知道它是如何来的。其实能带是一种近似的结果（可以看成一种近似），是周期边界条件（bloch函数）下的一种近似。先来看看一个最简单的问题，非周期体系有没有能带结构？答案是没有的，大家可以试试：
①建一个周期的晶胞②选择build菜单下的symmetry子菜单下的none periodic superstructure去掉周期边界条件性③看看还能够运行吗？运行（run）按钮变灰了，不能提交作业了。这说明什么问题？这说明这个CASTEP这个模块不能计算非周期的体系，另外可以参考MS中的DMOL模块，它可以计算非周期系统，虽然可以计算周期系统，但是仍不能计算能带，大家可以试试，看看property中的band structure能不能选上，一定不能！！^_^
从这里，我们可以得到一个结论，对于单个原子（分子、单胞）如果不加上周期边界条件，是无法获得能带结构的。所以计算小分子体系，或者采用团簇模型的朋友，这部分内容或许对你们没有帮助！那么，
非周期体系的态密度能够计算吗？这应该是能够计算的，曾经开到过文献采用团簇模型，计算出态密度的（phys. Rev. 上的文章）。
那么非周期体系为什么没有能带结构呢？
看一个例子：一个H2分子有能带吗？没有，因为它没有周期边界条件，也就是说在x，y，z方向上没有重复，所以它没有能带结构。那H2分子有什么东西呢？有两个轨道，两个
1s原子轨道，或者说两个轨道能级，它们成键参考右图。
再看另外一个例子：一维无限H原子链
H H H H H H
在一维无限H原子链体系中，产生了能带。
为什么在一维无限H原子链体系中能够产生能带呢？
因为，每一个H 原子有一个1s 轨道，由于在X 轴方向（H 原子周期排列的方向）引入周期边界条件，所以这个体系有无数（阿佛加得罗）个H1s 的轨道能级，这些具有相同能量的能级轨道处于简并的状态。如果两个相邻的H 原子之间距离较大，不能够成键，那么这无数个简并的能级将排成一条水平的直线，这条直线很长，无法画下来，那么我们只有压缩它，将他压缩到一个区间（[0,a π
]），这样每一个能级用一个点表示，由于点较多，看起
来好像形成了一条线，这样能带就形成了。
如果用函数的语言来描述，周期排列我们采用bloch 函数表示，我们解这些函数，就得到了一些k 矢量，在一维体系中k 矢量表示平移操作，k 的取值见下图，由于H1s 轨道能级有无数个，所以k 平移操作矢量就有无数个（注意k 是量
子化的），所以将它们压缩到[0,a π
]这个区间，就成了能带
结构中的横坐标，另外这个矢量也可以指向-X 方向，所以在[-
a π,0]这个区间能带的图像和[0,a π
]对称。
当H 原子之间的距离逐渐接近，它们的原子轨道要进行组合，形成一个成键分子轨道和一个反键分子轨道，那么原来能带是一条水平的直线，现在就要开始发生弯曲了（两个分子轨道能量不一样，导致能带发生弯曲），所以[0,a π
]这个区间，能带开始有带宽（散度），
随着H 原子的距离逐渐接近，可以预料，成键分子轨道和反键分子轨道的分裂越大，能带的带宽（最高能级－最低能级）越大。所以，相邻轨道之间的重叠越大，成键程度越大，带宽就越大。
另外值得一提的是，k 矢量还可以表示节点的数目。当k ＝0的时候表示什么呢？表示节点数为0，没有节点，所以k ＝0表示的H1s 轨道组合（组成分子轨道的原子轨道都带＋号）具有最大的成键，能量最低。随着k 值的增大，节点数逐渐增多，体系的能量上升，最后k ＝a π
时，H1s 轨道组合成的分子轨道能量最高（原子轨道为＋－＋－……交替）。所以H 原子链的能带结构是一条向上弯曲的曲线（能带），k ＝0能量最低，k ＝a π
时能量最高。
这里要特别注意，并不是k ＝0的时候节点数都为0，比如（＋－，＋－，＋－，……）这样的p 轨道，如果它们沿着X 轴周期排列，那么k ＝0的时候将具有最大的节点数，这时候形成的分子轨道将是能量最大的，随着k 的逐渐增加，节点数逐渐减小，所以这时候能带将向下伸展，这与H 原子链的情况刚好相反。
大家都知道，这样的H 原子链是不可能稳定的，最后都要变成H2分子，能带要消失，这是一个什么样的过程呢？在这个一维周期体系中轨道能级的数目假设为N （无限大），那么这个体系的电子数是多少呢？答案是N ，那么这些电子在这个能带上是如何分布的呢？当然是按照能量从低到高的顺序来填充的了，这样由于每一个轨道能够容纳2个电子，而这个能带只有较低能量的部分被填充（能带半充满），所以这时候要产生畸变（Peierls 畸变，即固体物理中的姜－泰勒效应），H 原子之间要产生相对振动（虚模，能量不稳定），以降低体系的能量，这样，H2分子就形成了，而能带也由于H2原子的形成破坏了周期条件，当这些H2分子不再沿着X 轴方向形成周期排列的时候，能带也就消失了（变成非周期体系）。
结论：一个原子的一个原子轨道在一维周期条件下将产生一条能带，能带的带宽取决于这些原子轨道的在周期方向上的成键强度，强度越大，带宽越大，成键越弱，带宽越小，如
果周期方向上没有成键，能带将是一条直线。另外能带是向上伸展还是向下伸展取决于原子轨道的特性，或者说是体系的拓补性质。
接下来我们看看，布里渊区里面的高对称点(G ，X ，F ，M 等)是怎么来的。
2、 布里渊区的高对称点
前面讲了一维周期条件下的布里渊区的能带是一条线，如果加上二维（X ，Y ）的周期边界条件，这些能带又会变成什么样呢？答案是一个面，由原来的线组成一个面。因为在一维情况下我们用一条能带来表示k 矢量（对称操作）和能级的关系，可以用E （k ）来表示，这构成第一布里渊区（即k 的取值范围[0,a π
]）。对于二维周期体系，我们需要两个平移矢
量k x 和k y ，所以能带可以用E （k x ，k y ）来表示，当k x ＝0时，变成E （0，k y ），得到一条能带（线，y 方向上与一维周期情况的能带类似）；当k y ＝0时，变成E （k x ，0），得到一条能带（线）。由于k x 和k y 是矢量，它们可以组合成另外一个矢
量，这个矢量不是沿着X 轴，也不是沿着Y 轴，实际上沿着该
矢量仍是能够得到一个能带（线）的，这样的矢量有很多，所
有的这些能带（线）将构成一个面。如果我们在做能带结构图
的时候，将能带结构按照二维的面画出来是很困难的，而三维
的情况更加困难，因为对称操作有很多，k 矢量的取值有很多，
所以一个可行的办法就是让k 的取值沿着一定的路径走，最后
回到起点。如右图（二维情况）。
这样，我们只要选择一些较高的对称点，就可以确定这个路径。比如二维的布里渊区是一个面，这个面上每一点与原点（G 点或Γ点）的连线都构成一个k 矢量，有一个k 矢量就有一个能级对应（E （k ））。所以，二维的能带结构是这个布里渊区上的一个平面（面积），如上图，按照Γ——>X ——>M ——>Γ这个路径走，就可以得到一个可以大致反映布里渊区上的能带平面的一个近似图，这就是二维的能带结构。具体的能带图的展开见下图。
三维的能带展开见下图：
前面讲了一维H 原子链的能带结构，提到这个体系具有一个能带，并且是向上弯曲的一个能带。这个能带是这样来的，H 原子链里面的基本单位（单体，单胞）是一个H 原子，每一个原子有一个原子轨道，即每一个基本单位有一个能级轨道，加上周期条件以后，这一系列轨道能就变为一个能带了。假如，现在以两个H 原子作为一个基本单位呢？能带结构又是如何的呢？这就是能带重叠的问题
3、 能带的折叠
如上图，将2个氢原子作为一个基本单位，这时候能带结构是什么样子的呢？很容易理解，根据前面的知识，单胞的原子轨道的数目决定了能带的数目，所以这样划分体系将有将有两个能带。
但是这个体系与前面1个H 原子周期链是一样的，只不过人为地进行了划分，能带结构就变了吗？是的，能带确实变了，那么能带将怎样变化呢？下图分别是1个H 原子为单胞和2个H 原子为单胞的能带结构。
第一个图可以清楚的看到，能带底部是成键的，能带的顶部是反键的，中间是非键的（成键与反键相当），能带向上伸展（弯曲）。在第二个图，原来的中部的非键轨道分别变成了两个能带的反键和成键
轨道（相同颜表示可以重叠成键，不同颜表示中间有一个节点）。实际上，这两个图是有关系的，能够反映相同的内容，首先，布里渊区从原来的[0,
a π]变为[0,a 2π
]，这也是可以理解的，原来的能带长度要变成原来的一半（因为周期方向上的单胞数减少一半，原来有N 个单胞，以2个H 原子为一个单胞后，单胞数变为一半，所以布里
渊区要减半）。其次，原来的能带在[0,a π
]展开，现在由于布里渊区减半，能带不能在[0,a 2π
]
这个小区间画出来，所以能带结构将产生折叠，由原来的一个能带变为2个能带，并且是以k ＝a 2π
为对称轴，将原来[a 2π，a π
]区间的能带折叠过去，所以就得到了2H 原子为单胞的
能带结构。如果我们将单胞取3个原子，或者取4个原子，体系的能带将如何变化？体系的能带分别变为
3条和4条能带，并且是2次或3次折叠原来的能带。
另外，上面的例子我们也可以从另外一个角度考虑问题，即2个H 原子组成一个H2分子，形成成键分子轨道和反键分子轨道，把这两个分子轨道看成一个“原子”的两个“原子轨道”，给这个“原子”加上周期条件以形成能带，这
样也得到2个能带。这种想法可以进行适当的扩充，
比如我们研究一个晶体，首先我们可以把单胞出来，
然后将单胞的轨道能级画出来，考虑它们之间的成键，
最后得到一组分子轨道，以这组分子轨道为基础，给
它们加上周期条件，形成能带。这样的想法什么时候
有用呢?在组合成分子轨道以后，这些分子轨道并不是
不变的，有可能和相邻晶胞的分子轨道再次组合，这
种情况在过渡金属氧化物中非常普遍，所以在考虑成
键的时候除了单胞内成键情况以外，还需要考虑晶胞间分子轨道的组合。
能带的折叠实际上是从不同的角度考虑问题而已，其实能带结构的本质还是一样的。如果我们理解了这点，那么我们在分析能带的时候，就不会笨到建立一个很大的超胞去研究其能带结构。所以我们要研究一个体系的能带结构时候，尽量选最小的单胞，以减少能带的数目。
不幸的是，有时候我们往往要计算掺杂的情况而不得不建立一个很大的单胞（超胞），而掺杂的比例又与超胞的大小有关。如果要做1/8掺杂，那么起码要将8个晶胞做成一个单胞（超胞），并将其中一个原子替换掉，从而实现掺杂。这样能带的数目将变为原来的8倍，这样我们计算出来的能带将是一团曲线交缠在一起，再也无法分开。
大家可以试试建立一个超胞，然后计算能带，看看能带结构如何！
3d走视图
所以能带结构的分析对于比较大的单胞（超胞）或者具有较多原子轨道的单胞将变得很麻烦。而且MS
中画出的能带结构很难看，如果我们要用origin 作图的画，又要反复的copy ＋paste 数据，所以建议在计算具有较多原子轨道（分子轨道）的体系时，尽量不用能带结构分析问题（虽然它很有用^_^）。用态密度和PDOS 就好多了。
另外，简单的提一下Fermi 能级的概念，有很多书都将最高占据分子轨道（HOMO ）的顶部作为Fermi 能级，所以我也采用这种观点，不知道大家是怎么理解的？
下面我简单讲讲我做的体系的能带结构，具体做法参加我提供的动画演示。
4、3维体系的能带——PbTiO3的能带结构
要分析能带，首先要知道每个能带是属于那个轨道的，所以首先要画出单胞的轨道能级图，下图是钛酸铅的轨道能级示意图，不是很准确（因为有一些轨道成键没有连线，仅给价带和导带部分连线了）。
在轨道能级图里面先将各个原子的轨道能级标出来（横轴的能量是参考了PDOS 加上的，没有也没有关系^_^）。接着考虑成键，主要是中心Ti 原子和O 原子成键（模型参考下图），有O 原子的2p 轨道组成价带，2s 轨道的能量最低Pb 的4s 轨道比O2p 轨道能量低，比O2s 轨道高。另外Ti 原子处于TiO6八面体中心，所以3d 能级要分裂成3下（3t 2g 态，3d xy ,3d xz ,3d
yz 单胞1单胞2
二次组合