电子结构计算DIY

电子结构计算DIY写这篇文章，主要是想解决大家心中类似这样的疑问：“我已经很熟悉第一原理的基本原理和公式了，但我不知道具体的电子结构计算过程是怎样通过电脑实现的，我输入了原子坐标，电脑是怎么求出电子结构呢？所有理论公式都是函数形式的，电脑是如何处理成数据的呢？”如果你也有同样的疑问，那么这篇文章将很适合你，我所要介绍的就是如何实现数值求解体系的电子结构，并通过一个简单的例子来展示如何手动计算电子结构。

但都只是科普性的介绍而已，对于长期使用软件却想了解其内部计算过程的同学会有些帮助，对于立志于编程的同志则有些过于初级了。

1．基本方法既然是科普性的介绍，那么我们就省去基本理论的介绍，简单的直接从体系本征方程开始吧 H E Ψ=Ψ（1）注意，在这里我们不想强调方程的具体形式，你可以把它看成KS 方程，也可以看成HF 方程，甚至看成薛定谔方程。

总之是体系波函数所满足的本征方程，其具体形式取决于所使用的理论方法，这里我们不想知道它是怎么来的，只想知道它是怎么解的。

无论如何，在我们的讨论中还是需要涉及一些理论和公式的。

首先我们看一下这个方程，体系波函数Ψ是完全未知的，直接让电脑描述这个函数是很不方便的，通常的做法是使用一组完备的基函数（f 1, f 2, f 3,…, fn ）把Ψ展开，1122n n c f c f c f Ψ=+++L（2）这里的基函数{f n }可以是平面波的组合，也可以是高斯型原子轨道函数，甚至可以像Dmol 软件那样：基函数不是解析而是数值形式的。

只要{f n }是完备的，就可以将体系的波函数展开成(2)的形式。

展开的好处是显而易见的：由于{f n }的形式是已知的，那么原本需要在整个空间内描述的Ψ，现在只需确定其展开系数{C n }就可以了。

具体的求解方法就是将(2)式代入本征方程(1)并求解。

《黑客帝国》中的matrix 是很玄的东西，与电脑紧密相关。

事实上的确如此，电脑处理矩阵问题的效率是很高的，我们这里的问题也很适合用矩阵解决，将(2)代入(1)后，本征方程可以写成下面的矩阵特征值方程的形式：1112111212222212n n n n nn n n H H H C C H H H C C E H H H C C     =  L L L L M M O M M M M M O M M M L L (3) 其中的矩阵元可以通过下式得到：ij j iH f H f = (4)其中的H 是体系的哈密顿算符。

对于确定的H 和f ，上式的积分利用电脑是可以方便的数值求解的（电脑处理积分的方法无非是密集取点、加和）。

上面提到了一定的理论推导，明显有悖于我科普的初衷，不过我们完全没有必要知道其中的理论细节，到目前只要记得：体系波函数可用基函数展开→可构造矩阵特征值方程→矩阵元可解，就可以了。

如果想对其中更深刻的原理，可以查看久期方程的知识。

需要强调的是(3)的解的意义，虽然其具体形式并不重要。

如果考研时突击的《线性代数》没有完全忘记的话，可以看出(3)就是n 维矩阵H 的特征方程，其有n 个特征值和相应的n 个特征向量。

这n 个特征值就是体系的电子能量分布情况，取其中最低的N 个能级是占据的。

而每个特征值E i 对应的特征向量{C i }（即[C i1, C i2, C i3,…, C in ]）呢？把它代入公式（2）就得到了本征值E i 所对应的电子波函数。

至此，我们已经可以利用电脑求解体系电子本征方程并得到电子能量以及波函数了。

但这仅仅是scf 过程的一步而已。

而scf 过程又是如何的呢？前面提到，对于确定的H 算符，可以得到矩阵元，进而得到特征方程，但实际上H 是未知的，是取决于体系电子状态的。

因此，真正的scf 过程是：构建初始波函数→ 利用波函数求解H 、构建矩阵方程并求解，得到新的波函数→ 比较新、旧波函数，若已经十分接近则scf 结束，否则重复上一步了解上述计算过程之后，我们首先可以搞清楚的一个事实就是第一原理计算量的问题。

可以看出，主要的计算出现在两个地方：1 积分计算矩阵元2 求解矩阵特征方程，而最主要的计算量是来自于第二点。

对于矩阵特征方程的求解过程，其计算量是O(n 3)的，即计算量与矩阵维度的3次方成正比，而矩阵的维度n 是由基函数的大小决定的。

因此，基函数是决定计算量的最重要因素。

具体到实际中，平面波法中的截断能、原子轨道法中的基组类型，都极大地影响着计算量。

（另外，体系的尺寸也影响着计算量。

原子轨道法的影响显而易见：原子数目直接决定着体系基组数目。

而平面波法的尺寸也影响着计算量，截断能只是决定了平面波基组中平面波波矢的取值范围，而波矢的取值密度是由晶胞大小决定的：晶胞越大，取值越密，基组也就越大。

因此晶胞尺寸影响着平面波法计算量。

这个可以通过比较原子轨道的dmol 和平面波的castep 处理团簇来看出：同样计算大超胞中的一个小团簇，可以发现dmol 的计算时间受超胞尺寸影响不大，而castep 的计算时间则随尺寸增加而疯狂增加）至于k 点对计算量的影响，很幸运的，k 点是可以“各个击破”的，就是说对于每一个k 点求解一次特征方程，而特征方程的维度是固定的，因此计算量随着k 点数目线性增加，这个影响远不如基组强烈。

这也就是为什么：（原胞+多k 点）的计算量远远小于（超胞+少k 点）的计算量。

好了，关于基本方法就说这么多，纸上得来终觉浅，下面我们通过小例子来自己手动算几个电子结构，以加深对上述过程的了解2．原子轨道线性组合（LCAO ）这一部分似乎涉及更多的理论公式，我会尽量“科普化”既然是手动计算，那么为了自己，啥也别说了，尽量简化问题吧！首先的简化就是我们将不进行scf 计算，只是一次性的求解本征方程，以此来感知一下电子结构的计算过程。

采用的进一步的简化方法是Slater 先生提出的LCAO 方法，具体理论可以参考 Phys. Rev. (94) 1498。

首先描述一下采用的基函数，我们采用原子轨道基函数。

对于非周期体系，设体系内位于R i 处的原子轨道为φn (r-R i )。

为了简化，我们用i 来区分不同原子，n 来标志不同的原子轨道,则基组可以简写成{φn,i }，即基组为所有原子、所有轨道的集合。

由于我们也想处理周期体系，作为简化，仅仅讨论一维的情况。

那么我们就要构建周期体系的基组，对于由N 个重复单元组成的一维周期结构，第j 个重复单元内的原子轨道可以写成φn(r-R i -R j )，其中R j 表示该单元的原点位置。

采用前面的简化方式，这个轨道可以表示为：φn,i (r-R j )根据Bloch 所做的工作，周期体系的Bloch sum 为：,,()()jik R n j n i j jk er R ψφ−•=−（5）可以对重复单元内每个原子的每个轨道都构建一个Bloch sum 。

对于周期体系，任一k 点的体系波函数均可以展开成Bloch sum{Ψn,i (k)}的线性组合（其中的理论根据，可以参考固体物理中的紧束缚方法以及布洛赫和部分的知识）这些Bloch sum 就是我们最终使用的基函数。

注意的是，非周期体系可以看成周期体系的特例，对于5式，若k=0，N=1,则基函数{Ψ}={φ}有了基组，只是解决了一个小问题，还有复杂的构建矩阵的过程等着我们呢。

为了方便讨论，我们用二个指标(n,j )来标识矩阵元，对某一k ，利用(4)可知矩阵元具有如下的形式：11222211(,)(,),,()()k n i n i n i n i H k H k ψψ=(6)由于基组中的e 指数项并不参与上式的积分，可以分离出来，因此所有矩阵元都含有2,21,1n i n i H φφ形式的积分。

而Slater 提出的LCAO 方法的核心就是：将这种典型的积分参数化，矩阵元通过这些参数的运算来得到，具体的细节我们将通过下面的实例来了解。

这里大概交待一下基本原理：比如我们已知原子i 1和i 2之间的积分2,21,1n i n i H φφ=A ，那么假设有另外一对原子i 3和i 4，其相对位置与i 1和i 2类似，那么积分4,43,3n i n i H φφ就可以利用A 值进行简单运算而得到。

因此对于一个复杂体系，只需一定数目的积分值，就可以将整个体系的哈密顿量得到。

（当然，这里涉及两中心近似，即积分2,21,1n i n i H φφ只与原子i 1和i 2的位置有关，与其他原子无关，这个假设在紧束缚近似下还是很不错的）上面就是我们所有的简化，又有很多公式推导，但同样地，我们完全没有必要知道其中的理论以及公式的具体细节，只要知道： 1 不进行scf 计算2 我们采用了原子轨道基组，是体系内所有原子以及原子轨道的组合（但我们根本不关心基组的具体形式）3 矩阵元不是通过计算积分得到，而是利用几个典型的积分值运算得到。

上面的这些简化使得我们的方法显得十分粗鄙，但出于手动计算的考虑，这些都是必须的，否则若需要计算多个积分才能实现，大家都不会感兴趣的。

另一方面，通过接下来的例子可以发现，这个方法还是比较靠谱的。

同时不要忘记，我们并不是用它来实际计算，只是用来了解电子结构计算的过程3．不多说了，开始实战了我要介绍的几个例子都是碳体系的，原因是碳体系电子结构比较简单，而且适合用LCAO 来描述。

在这个LCAO 模型中，我们仅仅考虑体系的π轨道电子相互作用，而忽略σ的作用。

因此这个模型只适合描述石墨体系。

另外一个假设是，我们只考虑最近邻原子的轨道相互作用，对于次近邻及以上的原子，认为之间的H 矩阵元是0。

首先把参数给出来，很简单的，只有一个数数：2.5。

这是石中最近邻原子的π轨道积分，也就是,2,1i i H ππφφ=2.5，其中i 1,i 2为近邻原子。

其它情况下积分都是零（注意，其实模型中还有另一个隐含的参数,1,1i i H ππφφ=0）。

我们从简单的非周期体系出发。

这种情况下基组就很简单了，体系每个原子只有一个π轨道πφ，体系基组就是{,n πφ} 首先处理一个如下图所示的体系，体系共含有22个C 原子，体系基组以及矩阵的维度n=22。

我们按照图中的方式对原子编号建立基组，并构建22维的哈密顿矩阵H ，根据公式(4)可知矩阵元,,ijj iH H ππφφ=，根据原子i 和j 是否近邻就可以判断矩阵元H ij 是2.5还是0。

例如，1号原子仅仅与2和8两个原子近邻，因此矩阵第一行中，仅第2和8列的矩阵元为2.5，其余均为0。

而2号原子只与1号和3 号原子近邻，因此第二行中，仅第1和3列的矩阵元为2.5，其余均为0，…按照这个方法，我们可以构建出整个矩阵来。

直到这里，我们才真正体会到了Slater 的LCAO 方法的优势：仅仅利用一个积分值就可以把整个矩阵解出来。