Dieict SCF的简单实现

zha23

在内存中存储N^4数目的双电子积分显然是不现实的，而用硬盘存储读取速度又是一大问题，现在的自洽场所用的基本都是直接方法。我自己写了一个简单的direct SCF，具体实现细节如下。1)积分的获取
我使用的c语言进行编程，电子积分库使用的为libcint，重叠积分和单电子算符积分在迭代开始之前计算并存储，双电子积分在迭代过程中计算。基组是从BSE上下载的。
2)线性代数库
我使用的为CBLAS和CLAPACK，这两个库由于都是用Fortran写的之后再用C写的接口，所以对于矩阵的排列些问题，C语言是行优先，而Fortran为列优先。
3)初猜
我写了一个extended Huckel进行初始猜想，但是实际的迭代效果并不理想，甚至比初始C=0的情况更糟。拓展Huckel方法通过原子的电子的电离能并使用最小基组计算H矩阵，对角化得到最小基组下的系数矩阵，但是现在大多使用的为劈裂和极化的基组，所以需要对得到的系数矩阵进行投影，实际就是对基进行变换，变换公式如下

其中P为最小基组和实际SCF计算基组之间的重叠积分，S为实际基组的重叠矩阵。
4)Fock矩阵计算
这一步是HF自洽场计算中最为费时的步骤，其时间复杂度为N^4。在直接SCF中每次都需要重新计算，那么太过于费时，所以可以用施瓦茨不等式提前判断积分大小的上限，同时direct SCF是在原来计算的fock矩阵上更新，所用不是密度矩阵而是前后两次迭代密度矩阵的差值。具体的判断公式如下：

计算完积分之后，需要对Fock矩阵进行更新，文献给出的公式如下：

我在这里费的时间最多，因为上面的公式其实并不适用所有的积分，其只在ijkl都不一样时才是正确的。计算双电子积分(ij|kl)，由对称性可以知道只需计算i>=j,k>=l,i+j>=kl的情况，实际计算量大约为所有双电子积分的1/8，所以很多积分不需要计算如(ji|kl)，但它们对Fock矩阵的贡献就需要在计算(ij|kl)时考虑，上述公式(4)中的系数4就是这么来的。对于库伦积分的情况较为简单，而交换积分的情况较为复杂，具体实现可以看代码，我的方法可能有些许问题，希望各位可以批评指正。
以上所有的密度矩阵都是针对RHF，D=C^{T}C，有的文章在讨论RHF时会使用D=2C^{T}C
5)能量计算
对Fock矩阵对角化可以得到一系列本征值e，由以下公式便可以计算RHF的能量

其时间复杂度为N^2
6)计算结果
对于乙烷使用6-31g*计算可以得到结果为-79.226225。

Freeman

其实只要利用施瓦茨不等式屏蔽掉绝对值很小的积分，利用双电子积分的八重对称性和分子的点群对称性避免重复计算等价的积分，存储成稀疏矩阵的形式，一般的计算节点其实是存得下中等大小的分子的双电子积分的。还不够的话就用密度拟合。
有的程序搞并行，可能把双电子积分复制了好几份分给每个进程，内存就不够用了，但是我猜只要把指向双电子积分所在的内存的指针分成好几份分给每个进程，就避免了复制双电子积分，也就节省了内存。

zha23

Freeman 发表于 2022-2-18 16:01
其实只要利用施瓦茨不等式屏蔽掉绝对值很小的积分，利用双电子积分的八重对称性和分子的点群对称性避免重复 ...

我没有进行并行化处理，我觉得不应该复制指针，而是直接复制数据，因为如果是多线程有访问冲突的问题，对于多进程其内存是不共享的。
我最大的问题是计算完积分后如何更新Fock矩阵，我计算的结果和gaussian的结果有些不同，并且用6-311g计算乙烷不收敛，而用6-31g**都可以。

Freeman

zha23 发表于 2022-2-18 16:14
我没有进行并行化处理，我觉得不应该复制指针，而是直接复制数据，因为如果是多线程有访问冲突的问题，对 ...

你说得对

zjxitcc

zha23 发表于 2022-2-18 16:14
我没有进行并行化处理，我觉得不应该复制指针，而是直接复制数据，因为如果是多线程有访问冲突的问题，对 ...

简单体系如果SCF不收敛，在2个能量值之间振荡的话，就要用DIIS了。