如何计算每一项双电子积分对Fock矩阵的贡献？

Freeman

Fock矩阵元写作：
F_ab=Hcore_ab+G_ab
其中
G_ab=Sigma(对c和d求和)P_cd[(ab|cd)-0.5(ad|cb)]
一个比较简单直接的求G_ab的方法就是四层循环，内二层循环是对于固定的(a,b)，使c和d分别独立地遍历所有原子轨道，以上式计算矩阵元G_ab，外二层循环是使a和b分别独立地遍历所有原子轨道，从而计算整个G矩阵。但是注意到(ab|cd)是有permutational symmetry的，也就是
(ab|cd)=(ab|dc)=(ba|cd)=(ba|dc)=(cd|ab)=(dc|ab)=(cd|ba)=(dc|ba)
所以为了节省计算资源，只要计算双电子积分总数的1/8即可。Libint给出的HF案例中，求G的方法是：
1、使用一套特殊的条件循环，不重样地算出每一个互相独立的双电子积分，即这些双电子积分互相不能用上面的恒等式连接；
2、计算每个双电子积分的多重度，比如说(11|11)只有它自己，多重度就是1，而(11|12)与(11|21)相等，所以多重度就是2；
3、计算每个不重样的双电子积分对G矩阵的贡献，也就相当于对F矩阵的贡献，具体写作：
G_12 += D_34 * (12|34) * 多重度
G_34 += D_12 * (12|34) * 多重度
G_13 -= 0.25 * D_24 * (12|34) * 多重度
G_24 -= 0.25 * D_13 * (12|34) * 多重度
G_14 -= 0.25 * D_23 * (12|34) * 多重度
G_23 -= 0.25 * D_14 * (12|34) * 多重度
其中D也是密度矩阵，和P的区别是，D是空间轨道的密度矩阵，而P是自旋轨道的密度矩阵。“+=”是C++的写法，表示累加，例如A+=B意为A=A+B。
我的问题是，上面6个G_矩阵元的公式是怎么来的呀？尤其是根据公式G_ab=Sigma(对c和d求和)P_cd[(ab|cd)-0.5(ad|cb)]，双电子积分对G的正贡献应该是负贡献的2倍，但是上面6个公式貌似体现出正贡献是负贡献的4倍了。这是咋回事儿？

以下是该HF示例文件，上面的一堆出现在第760行附近。 hartree-fock.cc (26.51 KB, 下载次数 Times of downloads: 10)

2021-8-7 23:30 上传 Uploaded

点击下载
Click to download

同样的问题我已经在Libint的issue页问过了，但是两天过去了，仍然没有获得答复。https://github.com/evaleev/libint/issues/214

wzkchem5

你可以把
G_12 += D_34 * (12|34) * 多重度
写作
G_12 += 0.5*(D_34+D_43) * (12|34) * 多重度
而对于交换项，比如以下面这个式子为例
G_13 -= 0.25 * D_24 * (12|34) * 多重度
如果把D_24换成D_42，(12|34)就要换成(14|32)了，这是一个不一样的积分，所以是另行计算的，所以交换项没办法像这样把D_24和D_42的贡献写在一起并把系数除以2。
这样就可以看出，库伦项的系数是交换项的2倍（0.5 vs. 0.25），不是4倍。

Freeman

wzkchem5 发表于 2021-8-8 01:24
你可以把
G_12 += D_34 * (12|34) * 多重度
写作

那么这个公式是怎么推出来的呢？
以G_12为例，G_12=Sigma(c,d)D_cd*[2*(12|cd)-(1d|c2)]。如果考虑(12|34)对它的贡献，(12|34)只对它有正贡献，具体来说就是以(12|34)和(12|43)的形式贡献在上式的(12|cd)项里，那么总共是4*D_34*(12|34)。但是如果用了G_12 += D_34 * (12|34) * 多重度 这个公式，因为(12|34)的多重度是8，所以G_12+= D_34*(12|34) *8。这样两者一个是4一个是8，不就矛盾了吗？

wzkchem5

Freeman 发表于 2021-8-8 04:27
那么这个公式是怎么推出来的呢？
以G_12为例，G_12=Sigma(c,d)D_cd*[2*(12|cd)-(1d|c2)]。如果考虑(12|3 ...

我估计代码最后还会有一个对称化，把(G_12+G_21)/2赋给G_12和G_21，这样就解释了系数差的那2倍

wangxubo

wzkchem5 发表于 2021-8-8 15:42
我估计代码最后还会有一个对称化，把(G_12+G_21)/2赋给G_12和G_21，这样就解释了系数差的那2倍

好估计，已经在源码和注释里看到了