开壳层Mayer键级计算求助

枫沨

最近闲来写了一个简单的Mayer键级的python程序，熟悉熟悉。不过对开壳层体系N2（正一价）我计算的结果跟multiwfn结果不一样。
奇怪的是，我计算出来的N2（+1）的键级竟然比中性的还要大，不清楚哪里做错了。

虽然分析是基于pyscf做的，但是我用mokit把pyscf轨道传给multiwfn计算的时候得到的结果也还是不一样。想在这里请教一下是什么原因导致的。
谢谢！


from pyscf import scf,gto,dft
import numpy as np
from mokit.lib import fchk

mol2 = gto.M(
    atom = '''
        N   0.00000000    0.00000000    0.54568900
        N   0.00000000    0.00000000   -0.54568900
            ''',
    unit = 'Ang',
    basis = 'ccpvtz',
    charge = 1,
    spin = 1,
)
# mol2.cart = False
# mol2.symmetry = False
mol2.build()

mf2 = dft.KS(mol2)
mf2.xc = 'b3lypg'
mf2.grids.atom_grid = (99,590)
mf2.newton()
mf2.kernel()

fchk(mf2, 'cation_N2_pyscf.fchk')

def get_aolist(mol):
    ao_label_list = mol.ao_labels()
    aos = np.zeros((mol.nao))
    for i in range(mol.nao):
        aos[i] = int(ao_label_list[i].split()[0])
    return aos

dm = mf2.make_rdm1()
ovlp = mol2.intor('int1e_ovlp')
aolist = get_aolist(mol2)

ps0 = np.matmul(dm[0], ovlp)
ps1 = np.matmul(dm[1], ovlp)

ii = np.where(aolist == 0)[0] # N1
jj = np.where(aolist == 1)[0] # N2

boma = 2 * np.einsum('ij,ji->', ps0[ii[:,None], jj], ps0[jj[:,None], ii])
bomb = 2 * np.einsum('ij,ji->', ps1[ii[:,None], jj], ps1[jj[:,None], ii])
print(boma, bomb, boma+bomb)
# 1.443662 1.723024 3.166686
# Multiwfn: 2.92540566

hebrewsnabla

你再仔细看看 Multiwfn的输出呢

1. 1
   
2. Calculating, please wait...
   
3. Bond orders with absolute value >=  0.050000
   
4. #    1:    1(N )    2(N ) Alpha:   1.443662 Beta:  1.723024 Total:  3.166686
   
5. 
   
6. Note: The "Total" bond orders shown above are more meaningful than the below ones. If you a
   
7. ut
   
8. 
   
9. Bond order from mixed alpha&beta density matrix >=  0.050000
   
10. #    1:         1(N )    2(N )    2.92540565
    
11. 
    
12. Total valences and free valences defined by Mayer:
    
13. Atom     1(N ) :    3.42915056    0.26246454
    
14. Atom     2(N ) :    3.42915056    0.26246454
    

复制代码

另外取ao的范围不用这么麻烦，可以这样

1. slice0, slice1 = mol2.aoslice_by_atom()
   
2. s0 = slice(slice0[2],slice0[3])
   
3. s1 = slice(slice1[2],slice1[3])
   
4. boma = 2 * np.einsum('ij,ji->',ps0[s0,s1], ps0[s1,s0])
   
5. bomb = 2 * np.einsum('ij,ji->',ps1[s0,s1], ps1[s1,s0])

复制代码

枫沨

hebrewsnabla 发表于 2024-9-6 17:51
你再仔细看看 Multiwfn的输出呢

啊，是我没仔细看输出信息。上面那个是alpha+beta的传统定义。下面那个mixed应该就是吧密度矩阵相加后用闭壳层方法算出来的结果。
还有多谢你的slice方法，确实比我的那个看起来整洁了很多！