matlab函数求和代码-Quantum-Chemistry:用于量子化学计算的代码集合,包括Hartree-Fock,耦合簇(

  • v3_168189
    了解作者
  • 55.8MB
    文件大小
  • zip
    文件格式
  • 0
    收藏次数
  • VIP专享
    资源类型
  • 0
    下载次数
  • 2022-04-06 04:30
    上传日期
matlab函数求和代码量子化学 基于Hartree-Fock和耦合簇(CC)方法论的用于量子化学计算的代码集合。 这些都是我自己为学习和测试目的而开发的代码。 存储库内容: CC_matlab 用Matlab编写的自旋积分耦合群集代码,实现了耦合群集(CC),完全重归一化(CR)CC,运动方程(EOM)CC和活动空间CC方法。 所有例程均经过矢量化处理,并且使用衍生自相似度转换的哈密顿量的中间体对方程进行了高度分解。 该代码旨在用作实现耦合集群方法的教学模板,以及用于测试新想法的沙箱。 使用从例如量子化学软件(例如GAMESS或Quantum Package 2.0)获得的一体和两体分子轨道积分,可以简单地运行代码。 它与RHF,UHF和ROHF参考状态兼容。 还包含SCF求解器,以便可以以自包含的方式生成分子轨道积分,尽管目前还没有对称适应性。 当前支持的计算包括: - RHF + RHF Analytical Gradients - CCSD + Left-CCSD - EOMCCSD + Left-EOMCCSD - CR-CC(2,3) + CR-EOMCC(2,3) - C
Quantum-Chemistry-master.zip
内容介绍
import numpy as np from solvers import diis import time import f90_updates #print(f90_updates.cc_loops.__doc__) def ccsdt(sys,ints,maxit=100,tol=1e-08,diis_size=6,shift=0.0): print('\n==================================++Entering CCSDT Routine++=================================\n') n1a = sys['Nocc_a'] * sys['Nunocc_a'] n1b = sys['Nocc_b'] * sys['Nunocc_b'] n2a = sys['Nocc_a'] ** 2 * sys['Nunocc_a'] ** 2 n2b = sys['Nocc_a'] * sys['Nocc_b'] * sys['Nunocc_a'] * sys['Nunocc_b'] n2c = sys['Nocc_b'] ** 2 * sys['Nunocc_b'] ** 2 n3a = sys['Nocc_a'] ** 3 * sys['Nunocc_a'] ** 3 n3b = sys['Nocc_a'] ** 2 * sys['Nocc_b'] * sys['Nunocc_a'] ** 2 * sys['Nunocc_b'] n3c = sys['Nocc_a'] * sys['Nocc_b'] ** 2 * sys['Nunocc_a'] * sys['Nunocc_b'] ** 2 n3d = sys['Nocc_b'] ** 3 * sys['Nunocc_b'] ** 3 ndim = n1a + n1b + n2a + n2b + n2c + n3a + n3b + n3c + n3d idx_1a = slice(0,n1a) idx_1b = slice(n1a,n1a+n1b) idx_2a = slice(n1a+n1b,n1a+n1b+n2a) idx_2b = slice(n1a+n1b+n2a,n1a+n1b+n2a+n2b) idx_2c = slice(n1a+n1b+n2a+n2b,n1a+n1b+n2a+n2b+n2c) idx_3a = slice(n1a+n1b+n2a+n2b+n2c,n1a+n1b+n2a+n2b+n2c+n3a) idx_3b = slice(n1a+n1b+n2a+n2b+n2c+n3a,n1a+n1b+n2a+n2b+n2c+n3a+n3b) idx_3c = slice(n1a+n1b+n2a+n2b+n2c+n3a+n3b,n1a+n1b+n2a+n2b+n2c+n3a+n3b+n3c) idx_3d = slice(n1a+n1b+n2a+n2b+n2c+n3a+n3b+n3c,n1a+n1b+n2a+n2b+n2c+n3a+n3b+n3c+n3d) cc_t = {} T = np.zeros(ndim) T_list = np.zeros((ndim,diis_size)) T_resid_list = np.zeros((ndim,diis_size)) T_old = np.zeros(ndim) # Jacobi/DIIS iterations it_micro = 0 flag_conv = False it_macro = 0 Ecorr_old = 0.0 t_start = time.time() print('Iteration Residuum deltaE Ecorr') print('=============================================================================') while it_micro < maxit: # store old T and get current diis dimensions T_old = T.copy() cc_t['t1a'] = np.reshape(T[idx_1a],(sys['Nunocc_a'],sys['Nocc_a'])) cc_t['t1b'] = np.reshape(T[idx_1b],(sys['Nunocc_b'],sys['Nocc_b'])) cc_t['t2a'] = np.reshape(T[idx_2a],(sys['Nunocc_a'],sys['Nunocc_a'],sys['Nocc_a'],sys['Nocc_a'])) cc_t['t2b'] = np.reshape(T[idx_2b],(sys['Nunocc_a'],sys['Nunocc_b'],sys['Nocc_a'],sys['Nocc_b'])) cc_t['t2c'] = np.reshape(T[idx_2c],(sys['Nunocc_b'],sys['Nunocc_b'],sys['Nocc_b'],sys['Nocc_b'])) cc_t['t3a'] = np.reshape(T[idx_3a],(sys['Nunocc_a'],sys['Nunocc_a'],sys['Nunocc_a'],sys['Nocc_a'],sys['Nocc_a'],sys['Nocc_a'])) cc_t['t3b'] = np.reshape(T[idx_3b],(sys['Nunocc_a'],sys['Nunocc_a'],sys['Nunocc_b'],sys['Nocc_a'],sys['Nocc_a'],sys['Nocc_b'])) cc_t['t3c'] = np.reshape(T[idx_3c],(sys['Nunocc_a'],sys['Nunocc_b'],sys['Nunocc_b'],sys['Nocc_a'],sys['Nocc_b'],sys['Nocc_b'])) cc_t['t3d'] = np.reshape(T[idx_3d],(sys['Nunocc_b'],sys['Nunocc_b'],sys['Nunocc_b'],sys['Nocc_b'],sys['Nocc_b'],sys['Nocc_b'])) # CC correlation energy Ecorr = calc_cc_energy(cc_t,ints) # update T1 cc_t = update_t1a(cc_t,ints,sys,shift) cc_t = update_t1b(cc_t,ints,sys,shift) # CCS intermediates H1A,H1B,H2A,H2B,H2C = get_ccs_intermediates(cc_t,ints,sys) # update T2 cc_t = update_t2a(cc_t,ints,H1A,H1B,H2A,H2B,H2C,sys,shift) cc_t = update_t2b(cc_t,ints,H1A,H1B,H2A,H2B,H2C,sys,shift) cc_t = update_t2c(cc_t,ints,H1A,H1B,H2A,H2B,H2C,sys,shift) # CCSD intermediates H1A,H1B,H2A,H2B,H2C = get_ccsd_intermediates(cc_t,ints,sys) # update T3 cc_t = update_t3a(cc_t,ints,H1A,H1B,H2A,H2B,H2C,sys,shift) cc_t = update_t3b(cc_t,ints,H1A,H1B,H2A,H2B,H2C,sys,shift) cc_t = update_t3c(cc_t,ints,H1A,H1B,H2A,H2B,H2C,sys,shift) cc_t = update_t3d(cc_t,ints,H1A,H1B,H2A,H2B,H2C,sys,shift) # store vectorized results T[idx_1a]= cc_t['t1a'].flatten() T[idx_1b] = cc_t['t1b'].flatten() T[idx_2a] = cc_t['t2a'].flatten() T[idx_2b] = cc_t['t2b'].flatten() T[idx_2c] = cc_t['t2c'].flatten() T[idx_3a] = cc_t['t3a'].flatten() T[idx_3b] = cc_t['t3b'].flatten() T[idx_3c] = cc_t['t3c'].flatten() T[idx_3d] = cc_t['t3d'].flatten() # build DIIS residual T_resid = T - T_old # change in Ecorr deltaE = Ecorr - Ecorr_old # check for exit condition resid = np.linalg.norm(T_resid) if resid < tol and abs(deltaE) < tol: flag_conv = True break # append trial and residual vectors to lists T_list[:,it_micro%diis_size] = T T_resid_list[:,it_micro%diis_size] = T_resid if it_micro%diis_size == 0 and it_micro > 1: it_macro = it_macro + 1 print('DIIS Cycle - {}'.format(it_macro)) T = diis(T_list,T_resid_list) print(' {} {:.10f} {:.10f} {:.10f}'.format(it_micro,resid,deltaE,Ecorr)) it_micro += 1 Ecorr_old = Ecorr t_end = time.time() minutes, seconds = divmod(t_end-t_start, 60) if flag_conv: print('CCSDT successfully converged! ({:0.2f}m {:0.2f}s)'.format(minutes,seconds)) print('') print('CCSDT Correlation Energy = {} Eh'.format(Ecorr)) print('CCSDT Total Energy = {} Eh'.format(Ecorr + ints['Escf'])) else: print('Failed to converge CCSDT in {} iterations'.format(maxit)) return cc_t, ints['Escf'] + Ecorr def update_t1a(cc_t,ints,sys,shift): vA = ints['vA'] vB = ints['vB'] vC = ints['vC'] fA = ints['fA'] fB = ints['fB'] t1a = cc_t['t1a'] t1b = cc_t['t1b'] t2a = cc_t['t2a'] t2b = cc_t['t2b'] t2c = cc_t['t2c'] t3a = cc_t['t3a'] t3b = cc_t['t3b'] t3c = cc_t['t3c'] t3d = cc_t['t3d'] chi1A_vv = 0.0 chi1A_vv += fA['vv'] chi1A_vv += np.einsum('anef,fn->ae',vA['vovv'],t1a,optimize=True) chi1A_vv += np.einsum('anef,fn->ae',vB['vovv'],t1b,optimize=True) chi1A_oo = 0.0 chi1A_oo += fA['oo'] chi1A_oo += np.einsum('mnif,fn->mi',vA['ooov'],t1a,optimize=True) chi1A_oo += np.einsum('mnif,fn->mi',vB['ooov'],t1b,optimize=True) h1A_ov = 0.0 h1A_ov += fA['ov'] h1A_ov += np.einsum('mnef,fn->me',vA['oovv'],t1a,optimize=True) h1A_ov += np.einsum('mnef,fn->me',vB['oovv'],t1b,optimize=True) h1B_ov = 0.0 h1B_ov += fB['ov'] h1B_ov += np.einsum('nmfe,fn->me',vB['oovv'],t1a,optimize=True) h1B_ov += np.einsum('mnef,fn->me',vC['oovv'],t1b,optimize=True) h1A_oo = 0.0 h1A_oo += chi1A_oo h1A_oo += np.einsum('me,ei->mi',h1A_ov,t1a,optimize=True) M11 = 0.0 M11 += fA['vo'] M11 -= np.einsum('mi,am->ai',h1A_oo,t1a,optimize=True) M11 += np.einsum('ae,ei->ai',chi1A_vv,t1a,optimize=True) M11 += np.einsum('anif,fn->ai',vA['voov'],t1a,optimize=True) M11 += np.einsum('anif,fn->ai',vB['voov'],t1b,optimize=True) h2A_ooov = vA['ooov'] + np.einsum('mnfe,fi->mnie',vA['oovv'],t1a,optimize=True) h2B_ooov = vB['ooov'] + np.einsum('mnfe,fi->mnie',vB['oovv'],t1a,optimize=True) h2A_vovv = vA['vovv'] - np.einsum('mnfe,an->amef',vA['oovv'],t1a,optimize=True) h2B_vovv = vB['vovv'] - np.einsum('nmef,an->amef',vB['oovv'],t1a,optimize=True) CCS_T2 = 0.0 CCS_T2 += np.einsum('me,aeim->ai',h1A_ov,t2a,optimize=True) CCS_T2 += np.einsum('me,aeim->ai',h1B_ov,t2b,optimize=True) CCS_T2 -= 0.5*np.einsum('mnif,afmn->ai',h2A_ooov,t2a,optimize=True) CCS_T2 -= np.einsum('mnif,afmn->ai',h2B_ooov,t2b,optimize=True) CCS_T2 += 0.5*np.einsum('anef,efin->ai',h2A_vovv,t2a,optimize=True) CCS_T2 += np.einsum('anef,efin->ai',h2B_vovv,t2b,optimize=True) X1A = M11 + CCS_T2 X1A += 0.25*np.einsum('mnef,aefimn->ai',vA['oovv'],t3a,optimize=True) X1A += np.einsum('mnef,aefimn->ai',vB['o
评论
    相关推荐