计算化学公社

标题: 自写NOCV的python脚本计算错误，求意见 [打印本页]
作者
Author: flyingchow    时间: 2023-6-27 05:32
标题: 自写NOCV的python脚本计算错误，求意见
本帖最后由 flyingchow 于 2023-6-27 05:37 编辑

在看了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)
复制代码
请老师帮忙看看，到底问题出在哪里。谢谢大家了。





欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3