PySCF中双电子积分索引问题：eri1形状与<ij||kl>积分提取

我正在开发一个改进版的FCI程序包，学习行列式索引、Slater规则应用等内容，目前在双电子积分环节遇到问题。我需要获取基组对应的积分集合，用于系统的Slater-Condon规则求解，已完成部分测试并成功得到单电子积分，但疑惑为何eri1的形状为(100,100)，以及如何从该二维数组中提取特定的<ij||kl>积分。

测试代码 #

import pyscf as ps
import pyscf.mcscf as mc
import pyscf.tools.molden as md
import pyscf.lo as lo
import scipy
import pyscf.scf.ghf as ghf
import pyscf.ao2mo as ao2mo
import numpy as np

print("yes")

mol_h2 = ps.gto.M(atom = 'H 0 1 0; H 0 0 1', basis = 'ccpvdz', symmetry=False)

uhf_h2 = ps.scf.UHF(mol_h2)
uhf_h2.kernel()
print(mol_h2.intor('int1e_kin').shape)
print(mol_h2.intor('int1e_nuc').shape)
print(mol_h2.intor('int1e_ovlp').shape)
print(mol_h2.intor('int2e').shape)
print(uhf_h2.mo_coeff.shape)
print(uhf_h2.mo_occ.shape)

eri1 = mol_h2.ao2mo(uhf_h2.mo_coeff[0], uhf_h2.mo_coeff[0], aosym = 's1', compact=False)
print(eri1.shape)

H=np.zeros([2,2])

N = 10
sum=0
for i in np.arange(0,N):
    if(uhf_h2.mo_occ[0,i]==1):
        for j in np.arange(0,N):
            for k in np.arange(0,N):
                #sum = sum+ uhf_h2.mo_coeff[0,i,j]*uhf_h2.mo_coeff[0,i,k]*(mol_h2.intor('int1e_kin')[j,k]+mol_h2.intor('int1e_nuc')[j,k])
                sum = sum+ uhf_h2.mo_coeff[0,j,i]*uhf_h2.mo_coeff[0,k,i]*(mol_h2.intor('int1e_kin')[j,k]+mol_h2.intor('int1e_nuc')[j,k])
    if(uhf_h2.mo_occ[1,i]==1):
        for j in np.arange(0,N):
            for k in np.arange(0,N):
                #sum = sum+ uhf_h2.mo_coeff[0,i,j]*uhf_h2.mo_coeff[0,i,k]*(mol_h2.intor('int1e_kin')[j,k]+mol_h2.intor('int1e_nuc')[j,k])
                sum = sum+ uhf_h2.mo_coeff[1,j,i]*uhf_h2.mo_coeff[1,k,i]*(mol_h2.intor('int1e_kin')[j,k]+mol_h2.intor('int1e_nuc')[j,k])

输出结果 #

yes
converged SCF energy = -1.01844561914196  <S^2> = 1.5343282e-13  2S+1 = 1
(10, 10)
(10, 10)
(10, 10)
(10, 10, 10, 10)
(2, 10, 10)
(2, 10)
(100, 100)

问题解答 #

一、eri1形状为(100,100)的原因 #

你的体系是两个H原子，使用cc-pVDZ基组后总共包含10个原子轨道（AO），对应10个分子轨道（MO）。调用ao2mo时，参数compact=False和aosym='s1'指定了不压缩积分、不利用对称性优化存储。此时pyscf会将四阶双电子积分（原本形状为(10,10,10,10)）的前两个指标(i,j)合并为一个一维索引，后两个指标(k,l)也合并为一个一维索引：

• 每个合并后的索引长度为10*10=100（i和j各有10种取值）
• 最终得到的二维数组对应(ij, kl)的存储格式，其中ij = i * nmo + j，kl = k * nmo + l（nmo=10为MO轨道数）

二、提取特定<ij||kl>积分的方法 #

Slater-Condon规则中用到的<ij||kl>是反对称化双电子积分，定义为：
$$<ij||kl> = <ij|kl> - <il|kj>$$

你可以通过两种方式获取：

1. 直接从(100,100)的eri1数组计算：
   先将i,j和k,l转换为合并后的一维索引，再代入公式计算：

   nmo = 10
   # 定义要提取的轨道索引i,j,k,l
   i, j, k, l = 0, 1, 2, 3
   idx_ij = i * nmo + j
   idx_kl = k * nmo + l
   idx_il = i * nmo + l
   idx_kj = k * nmo + j
   ijkl_antisym = eri1[idx_ij, idx_kl] - eri1[idx_il, idx_kj]
   

2. 将eri1恢复为四阶张量后直接索引：
   使用pyscf的ao2mo.restore函数把二维数组转回四阶张量，操作更直观：

   from pyscf import ao2mo
   nmo = 10
   eri_4d = ao2mo.restore(1, eri1, nmo)  # 参数1表示不利用对称性，恢复为(10,10,10,10)的四阶数组
   # 定义要提取的轨道索引i,j,k,l
   i, j, k, l = 0, 1, 2, 3
   ijkl_antisym = eri_4d[i,j,k,l] - eri_4d[i,l,k,j]
   

内容的提问来源于stack exchange，提问作者Elijah123