自写NOCV的python脚本计算错误，求意见

flyingchow

在看了Multiwfn的NOCV源码之后，我就想用python计算NOCV的eigenvalue，计算过程和Multiwfn是一样的。但是得不到最后成对的eigenvalue。
我现在把这个错误的脚本发到论坛上来，各位老师检查一下，看看哪里能是我搞错了。

我原计划是支持pyscf和turbomole两种软件。
构建CObasapro

1. # build diagonal block coefficient matrix for fragments A and B
   
2.     def get_coeff_pro(self):
   
3.         # obtain the coefficient matrix of fragments A and B
   
4.         fragA_coeff = self.mol1_coeff
   
5.         fragB_coeff = self.mol2_coeff
   
6.         zero_like_A = np.zeros((len(fragA_coeff),len(fragB_coeff)))
   
7.         zero_like_B = np.zeros((len(fragB_coeff),len(fragA_coeff)))
   
8.         # build the full coefficient matrix for fragments A and B
   
9.         self.coeff_pro = np.block([[fragA_coeff, zero_like_A],
   
10.                                 [zero_like_B, fragB_coeff]])

复制代码

计算nocc/norb 总轨道和occupied轨道 以及occupied轨道在CObasapro中的列数

1. def get_occupied(self):
   
2.         # obtain the column number of occupied orbitals in s_frag and coeff_frz
   
3.         self.nocc = self.mol_3.nelectron//2
   
4.         norb = self.mol_3.nao_nr()
   
5.         nocc_A = self.mol_1.nelectron//2
   
6.         norb_A = self.mol_1.nao_nr()
   
7.         nocc_B = self.mol_2.nelectron//2
   
8.         norb_B = self.mol_2.nao_nr()
   
9.         self.occupied_rows = []
   
10.         if norb == norb_A + norb_B:
    
11.             for i in range(nocc_A):
    
12.                 self.occupied_rows.append(i)
    
13.             # starting from occ_A, ending at occ_A+nocc_B
    
14.             for i in range(norb_A,norb_A+nocc_B):
    
15.                 self.occupied_rows.append(i)
    
16.             print (f'{norb} basis sets were used and {self.nocc} orbitals are occupied')
    
17.             print(f'{norb_A} basis sets were used for fragment A and {nocc_A} orbitals are occupied')
    
18.             print(f'{norb_B} basis sets were used for fragment B and {nocc_B} orbitals are occupied')

复制代码

计算Cobasa_frz和P_diff

1.     def get_coeff_frz(self):
   
2.         # obtain the overlap matrix of mol3 (AB combined)
   
3.         self.s = self.mol_3.intor_symmetric('int1e_ovlp')
   
4.         # calculate MO overlap matrix from AO overlap matrix
   
5.         self.s_frag_total = matmul(self.coeff_pro.T,matmul(self.s,self.coeff_pro))
   
6.         self.s_frag = self.s_frag_total[self.occupied_rows, :][:, self.occupied_rows]
   
7.         # get eigenvalues and eigenvectors of s_frag
   
8.         eigval_s_frag,eigvec_s_frag = np.linalg.eigh(self.s_frag)
   
9.         # calculate eigenvalues inverse square root and get S^(-1/2)
   
10.         eigval_inv_sqrt = np.diag(1 / np.sqrt(eigval_s_frag))
    
11.         self.s_frag_inv_sqrt = matmul(eigvec_s_frag,matmul(eigval_inv_sqrt,eigvec_s_frag.T))
    
12.         # get coefficient matrix of frozen orbitals
    
13.         self.coeff_frz = np.einsum('ij,jk->ik',self.coeff_pro[:,self.occupied_rows],self.s_frag_inv_sqrt)
    
14.         # get density matrix of frozen orbitals
    
15.         self.den_frz = np.einsum('pi,qi->pq', self.coeff_frz, np.conj(self.coeff_frz))*2
    
16.     def get_nocv_eig(self):
    
17.        self.den_deform = self.den_total - self.den_frz
    
18.        self.den_deform = np.matmul(self.den_deform, self.s)
    
19.        self.nocv_eigval,self.nocv_eigvec = np.linalg.eigh(self.den_deform)

复制代码

请老师帮忙看看，到底问题出在哪里。谢谢大家了。