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

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