求助一下libcint库算双电子积分和pyscf给出结果不一致的问题

海棠

这个是我自己编写低效率Obara-Saika方案双电子积分的结果：

i=  1 j=  1 k=  1 l=  1 int=  0.165539885177D+01
i=  2 j=  1 k=  1 l=  1 int=  0.105390230958D+00
i=  2 j=  1 k=  2 l=  1 int=  0.130812876581D-01
i=  2 j=  2 k=  1 l=  1 int=  0.370298294556D+00

这个是我用libcint库得到的双电子积分结果：

I=  1 J=  1 K=  1 L=  1 Int=  0.165046544094D+01
I=  2 J=  1 K=  1 L=  1 Int=  0.105422335371D+00
I=  2 J=  1 K=  2 L=  1 Int=  0.130757373222D-01
I=  2 J=  2 K=  1 L=  1 Int=  0.370292716483D+00

这个是Gaussian09的输入文件：

#p hf/gen iop(3/33=2) nosymm IOp(3/33=3) IOp(5/33=4) extralinks(L316)
# scf=conventional noraff

Li2 spherical harmonic basis

0 1
Li                0.0   0.0    0.00000000
Li                0.0   0.0    2.7252628323500515

Li     0
S    6   1.00
      0.6424189150D+03       0.2142607810D-02
      0.9679851530D+02       0.1620887150D-01
      0.2209112120D+02       0.7731557250D-01
      0.6201070250D+01       0.2457860520D+00
      0.1935117680D+01       0.4701890040D+00
      0.6367357890D+00       0.3454708450D+00
SP   3   1.00
      0.2324918408D+01      -0.3509174574D-01       0.8941508043D-02
      0.6324303556D+00      -0.1912328431D+00       0.1410094640D+00
      0.7905343475D-01       0.1083987795D+01       0.9453636953D+00
SP   1   1.00
      0.3596197175D-01       0.1000000000D+01       0.1000000000D+01
SP   1   1.00
      0.7400000000D-02       0.1000000000D+01       0.1000000000D+01
D    1   1.00
      0.2000000000D+00       1.0000000
****

这个是Gaussian09给出的双电子积分结果：

I=  2 J=  1 K=  2 L=  1 Int=  0.130812876572D-01
I=  2 J=  2 K=  1 L=  1 Int=  0.370298294532D+00
I=  2 J=  1 K=  1 L=  1 Int=  0.105390230966D+00
I=  1 J=  1 K=  1 L=  1 Int=  0.165539885212D+01

这个是pyscf给出的结果：

>>> from pyscf import gto
>>> mol_li2=gto.M(atom='Li 0 0 0;Li 0 0 2.7252628323500515',basis='6-31+g*')
>>> eri_4fold_ao=mol_li2.intor('int2e_sph')
>>> eri_4fold_ao[0,0,0,0]
1.655398912142564
>>> eri_4fold_ao[1,0,0,0]
0.10539019493053174
>>> eri_4fold_ao[1,0,1,0]
0.013081278619072151
>>> eri_4fold_ao[1,1,0,0]
0.3702982386189271

可以看到pyscf给出的双电子积分和Gaussian09给出的，比较相近（自己写的低效率程序也可以与二者很好地接近）。但我自己写的程序调用libcint库却得到了差异比较大的结果。

另，在调用libcint前给数组env赋值，如果按libcint/example里的示例程序，先给出原子坐标，再给出基的指数和收缩系数，那么得到如上的结果；如果env里先给出基的指数和收缩系数，再给出原子坐标，则双电子积分差别更大。由于那种情况的结果没有保存，就不附在此处了。这些信息在env数组的位置应当由数组atm和bas指定，它们哪个在前面哪个在后面本应该没有影响才对啊？

所以我求助一下，我使用libcint的时候出了这么多岔子，可能的错误出现在哪里？

为了说明我自认为基组信息没有读错，附自写程序用libcint库算出的Core Hamiltonian矩阵与Gaussian09的对比（部分）。

Gaussian09：

      1 -0.502779D+01
      2 -0.613493D+00 -0.154692D+01
      3  0.000000D+00  0.000000D+00 -0.129947D+01
      4  0.000000D+00  0.000000D+00  0.000000D+00 -0.129947D+01
      5 -0.148923D-01 -0.201703D+00  0.000000D+00  0.000000D+00 -0.146475D+01
      6 -0.842818D+00 -0.138291D+01  0.000000D+00  0.000000D+00 -0.214278D+00
      7  0.000000D+00  0.000000D+00 -0.929379D+00  0.000000D+00  0.000000D+00
      8  0.000000D+00  0.000000D+00  0.000000D+00 -0.929379D+00  0.000000D+00
      9 -0.408435D-02 -0.159307D+00  0.000000D+00  0.000000D+00 -0.108801D+01
     10 -0.269105D+00 -0.628740D+00  0.000000D+00  0.000000D+00 -0.117735D+00
     11  0.000000D+00  0.000000D+00 -0.248969D+00  0.000000D+00  0.000000D+00
     12  0.000000D+00  0.000000D+00  0.000000D+00 -0.248969D+00  0.000000D+00
     13 -0.622566D-03 -0.403141D-01  0.000000D+00  0.000000D+00 -0.293966D+00
     14 -0.331608D-02 -0.879144D-01  0.000000D+00  0.000000D+00 -0.203405D+00
     15  0.000000D+00  0.000000D+00 -0.110811D+00  0.000000D+00  0.000000D+00
     16  0.000000D+00  0.000000D+00  0.000000D+00 -0.110811D+00  0.000000D+00
     17  0.000000D+00  0.000000D+00  0.000000D+00  0.000000D+00  0.000000D+00
     18  0.000000D+00  0.000000D+00  0.000000D+00  0.000000D+00  0.000000D+00
     19 -0.184754D-03 -0.225257D+00  0.000000D+00  0.000000D+00 -0.523731D+00
     20 -0.225257D+00 -0.683483D+00  0.000000D+00  0.000000D+00 -0.867866D+00
     21  0.000000D+00  0.000000D+00 -0.461770D+00  0.000000D+00  0.000000D+00
     22  0.000000D+00  0.000000D+00  0.000000D+00 -0.461770D+00  0.000000D+00
     23  0.523731D+00  0.867866D+00  0.000000D+00  0.000000D+00  0.663913D+00
     24 -0.337799D+00 -0.803886D+00  0.000000D+00  0.000000D+00 -0.728719D+00
     25  0.000000D+00  0.000000D+00 -0.517956D+00  0.000000D+00  0.000000D+00
     26  0.000000D+00  0.000000D+00  0.000000D+00 -0.517956D+00  0.000000D+00
     27  0.633924D+00  0.984422D+00  0.000000D+00  0.000000D+00  0.282532D+00
     28 -0.221720D+00 -0.540657D+00  0.000000D+00  0.000000D+00 -0.210284D+00
     29  0.000000D+00  0.000000D+00 -0.212240D+00  0.000000D+00  0.000000D+00
     30  0.000000D+00  0.000000D+00  0.000000D+00 -0.212240D+00  0.000000D+00
     31  0.194385D+00  0.403085D+00  0.000000D+00  0.000000D+00 -0.890999D-01
     32 -0.251864D+00 -0.374486D+00  0.000000D+00  0.000000D+00 -0.315558D+00
     33  0.000000D+00  0.000000D+00  0.339238D+00  0.000000D+00  0.000000D+00
     34  0.000000D+00  0.000000D+00  0.000000D+00  0.339238D+00  0.000000D+00
     35  0.000000D+00  0.000000D+00  0.000000D+00  0.000000D+00  0.000000D+00
     36  0.000000D+00  0.000000D+00  0.000000D+00  0.000000D+00  0.000000D+00

海棠

另外，我使用libcint库只能是编译出libcint.so libcint.so.5 libcint.so.5.1.3三个文件后自行移动到/usr/lib/里才能使用，使用ifort -L命令没有效果。不知道上述错误是否和我没有正常链接动态库有关？

hebrewsnabla

这些信息在env数组的位置应当由数组atm和bas指定，它们哪个在前面哪个在后面本应该没有影响才对啊？

并不是这样，至少env的前二十个数是不能随便动的，见这里sunqm的讨论 https://github.com/sunqm/libcint/issues/70 。

不过，我觉得并不一定是这个问题导致的。pyscf生成env的时候实际上对基组做了归一化，见 https://github.com/pyscf/pyscf/blob/master/pyscf/gto/mole.py 的 make_bas_env 函数。如果你没有做，那结果自然和pyscf的不同。

如果你想和pyscf一致，那必须仔细阅读pyscf生成env的代码。但是，不一致也没有关系，你这样写也是能用的。

hebrewsnabla

海棠 发表于 2022-7-2 13:34
另外，我使用libcint库只能是编译出libcint.so libcint.so.5 libcint.so.5.1.3三个文件后自行移动到/usr/li ...

无关

海棠

hebrewsnabla 发表于 2022-7-2 14:16
并不是这样，至少env的前二十个数是不能随便动的，见这里sunqm的讨论 https://github.com/sunqm/libcint/ ...

非常感谢！把前20个元素空出来之后得到了合理结果。