计算化学公社

标题: 开壳层Mayer键级计算求助 [打印本页]

作者
Author: 枫沨 时间: 2024-9-6 16:49
标题: 开壳层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

作者
Author: hebrewsnabla 时间: 2024-9-6 17:51
你再仔细看看 Multiwfn的输出呢

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

复制代码

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

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

复制代码

作者
Author: 枫沨 时间: 2024-9-6 19:06

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

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

欢迎光临计算化学公社 (http://bbs.keinsci.com/)