From b6c74723bccded384900da2a6109d17ec0930701 Mon Sep 17 00:00:00 2001 From: Unknown Date: Tue, 24 Jul 2018 16:10:00 -0700 Subject: [PATCH 1/7] Cosmetics, unused import --- greens_function.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/greens_function.py b/greens_function.py index c8b0317..f64bcce 100644 --- a/greens_function.py +++ b/greens_function.py @@ -1,4 +1,3 @@ -import collections import numpy as np import scipy import gminres @@ -6,10 +5,12 @@ import pyscf import pyscf.cc from pyscf.cc.eom_rccsd import amplitudes_to_vector_ip, amplitudes_to_vector_ea + ################### # EA Greens # ################### + def greens_b_vector_ea_rhf(cc,p): nocc, nvir = cc.t1.shape ds_type = cc.t1.dtype @@ -23,6 +24,7 @@ def greens_b_vector_ea_rhf(cc,p): vector1[ p-nocc ] = 1.0 return amplitudes_to_vector_ea(vector1,vector2) + def greens_e_vector_ea_rhf(cc,p): nocc, nvir = cc.t1.shape ds_type = cc.t1.dtype @@ -49,6 +51,7 @@ def greens_e_vector_ea_rhf(cc,p): # IP Greens # ################### + def greens_b_vector_ip_rhf(cc,p): nocc, nvir = cc.t1.shape vector1 = np.zeros((nocc),dtype=complex) @@ -60,6 +63,7 @@ def greens_b_vector_ip_rhf(cc,p): vector2 += cc.t2[:,:,:,p-nocc] return amplitudes_to_vector_ip(vector1,vector2) + def greens_e_vector_ip_rhf(cc,p): nocc, nvir = cc.t1.shape vector1 = np.zeros((nocc),dtype=complex) @@ -80,15 +84,18 @@ def greens_e_vector_ip_rhf(cc,p): return amplitudes_to_vector_ip(vector1,vector2) + def greens_func_multiply(ham,vector,linear_part,args=None): return np.array(ham(vector) + (linear_part)*vector) + def initial_ip_guess(cc): nocc, nvir = cc.t1.shape vector1 = np.zeros((nocc),dtype=complex) vector2 = np.zeros((nocc,nocc,nvir),dtype=complex) return amplitudes_to_vector_ip(vector1,vector2) + def initial_ea_guess(cc): nocc, nvir = cc.t1.shape vector1 = np.zeros((nvir),dtype=complex) @@ -100,15 +107,16 @@ def ip_shape(cc): nocc, nvir = cc.t1.shape return nocc + nocc*nocc*nvir + def ea_shape(cc): nocc, nvir = cc.t1.shape return nvir + nocc*nvir*nvir + class greens_function: def __init__(self, verbose=0): self.verbose = verbose - def td_ip_ao(self,cc,ps,times,mo_coeff,re_im="re",tol=1.e-5): """ E0: total CC gs energy From edc0acbfa0b071085d1fe5a12038ea61586c471b Mon Sep 17 00:00:00 2001 From: Unknown Date: Tue, 24 Jul 2018 17:53:21 -0700 Subject: [PATCH 2/7] Reformatted code for readablity (pycharm) --- greens_function.py | 393 +++++++++++++++++++++++---------------------- 1 file changed, 201 insertions(+), 192 deletions(-) diff --git a/greens_function.py b/greens_function.py index f64bcce..d37cce5 100644 --- a/greens_function.py +++ b/greens_function.py @@ -6,118 +6,120 @@ import pyscf.cc from pyscf.cc.eom_rccsd import amplitudes_to_vector_ip, amplitudes_to_vector_ea + ################### # EA Greens # ################### -def greens_b_vector_ea_rhf(cc,p): +def greens_b_vector_ea_rhf(cc, p): nocc, nvir = cc.t1.shape ds_type = cc.t1.dtype - vector1 = np.zeros((nvir),dtype=ds_type) - vector2 = np.zeros((nocc,nvir,nvir),dtype=ds_type) + vector1 = np.zeros((nvir), dtype=ds_type) + vector2 = np.zeros((nocc, nvir, nvir), dtype=ds_type) if p < nocc: # Changed both to minus - vector1 += -cc.t1[p,:] - vector2 += -cc.t2[p,:,:,:] + vector1 += -cc.t1[p, :] + vector2 += -cc.t2[p, :, :, :] else: - vector1[ p-nocc ] = 1.0 - return amplitudes_to_vector_ea(vector1,vector2) + vector1[p - nocc] = 1.0 + return amplitudes_to_vector_ea(vector1, vector2) -def greens_e_vector_ea_rhf(cc,p): +def greens_e_vector_ea_rhf(cc, p): nocc, nvir = cc.t1.shape ds_type = cc.t1.dtype - vector1 = np.zeros((nvir),dtype=ds_type) - vector2 = np.zeros((nocc,nvir,nvir),dtype=ds_type) + vector1 = np.zeros((nvir), dtype=ds_type) + vector2 = np.zeros((nocc, nvir, nvir), dtype=ds_type) if p < nocc: # Changed both to plus - vector1 += cc.l1[p,:] - vector2 += (2*cc.l2[p,:,:,:] - cc.l2[:,p,:,:]) + vector1 += cc.l1[p, :] + vector2 += (2 * cc.l2[p, :, :, :] - cc.l2[:, p, :, :]) pass else: - vector1[ p-nocc ] = -1.0 - vector1 += np.einsum('ia,i->a', cc.l1, cc.t1[:,p-nocc]) - vector1 += 2*np.einsum('klca,klc->a', cc.l2, cc.t2[:,:,:,p-nocc]) - vector1 -= np.einsum('klca,lkc->a', cc.l2, cc.t2[:,:,:,p-nocc]) + vector1[p - nocc] = -1.0 + vector1 += np.einsum('ia,i->a', cc.l1, cc.t1[:, p - nocc]) + vector1 += 2 * np.einsum('klca,klc->a', cc.l2, cc.t2[:, :, :, p - nocc]) + vector1 -= np.einsum('klca,lkc->a', cc.l2, cc.t2[:, :, :, p - nocc]) + + vector2[:, p - nocc, :] += -2. * cc.l1 + vector2[:, :, p - nocc] += cc.l1 + vector2 += 2 * np.einsum('k,jkba->jab', cc.t1[:, p - nocc], cc.l2) + vector2 -= np.einsum('k,jkab->jab', cc.t1[:, p - nocc], cc.l2) + return amplitudes_to_vector_ea(vector1, vector2) - vector2[:,p-nocc,:] += -2.*cc.l1 - vector2[:,:,p-nocc] += cc.l1 - vector2 += 2*np.einsum('k,jkba->jab', cc.t1[:,p-nocc], cc.l2) - vector2 -= np.einsum('k,jkab->jab', cc.t1[:,p-nocc], cc.l2) - return amplitudes_to_vector_ea(vector1,vector2) ################### # IP Greens # ################### -def greens_b_vector_ip_rhf(cc,p): +def greens_b_vector_ip_rhf(cc, p): nocc, nvir = cc.t1.shape - vector1 = np.zeros((nocc),dtype=complex) - vector2 = np.zeros((nocc,nocc,nvir),dtype=complex) + vector1 = np.zeros((nocc), dtype=complex) + vector2 = np.zeros((nocc, nocc, nvir), dtype=complex) if p < nocc: - vector1[ p ] = 1.0 + vector1[p] = 1.0 else: - vector1 += cc.t1[:,p-nocc] - vector2 += cc.t2[:,:,:,p-nocc] - return amplitudes_to_vector_ip(vector1,vector2) + vector1 += cc.t1[:, p - nocc] + vector2 += cc.t2[:, :, :, p - nocc] + return amplitudes_to_vector_ip(vector1, vector2) -def greens_e_vector_ip_rhf(cc,p): +def greens_e_vector_ip_rhf(cc, p): nocc, nvir = cc.t1.shape - vector1 = np.zeros((nocc),dtype=complex) - vector2 = np.zeros((nocc,nocc,nvir),dtype=complex) + vector1 = np.zeros((nocc), dtype=complex) + vector2 = np.zeros((nocc, nocc, nvir), dtype=complex) if p < nocc: - vector1[ p ] = -1.0 - vector1 += np.einsum('ia,a->i', cc.l1, cc.t1[p,:]) - vector1 += 2*np.einsum('ilcd,lcd->i', cc.l2, cc.t2[p,:,:,:]) - vector1 -= np.einsum('ilcd,ldc->i', cc.l2, cc.t2[p,:,:,:]) - - vector2[p,:,:] += -2.*cc.l1 - vector2[:,p,:] += cc.l1 - vector2 += 2*np.einsum('c,ijcb->ijb', cc.t1[p,:], cc.l2) - vector2 -= np.einsum('c,jicb->ijb', cc.t1[p,:], cc.l2) + vector1[p] = -1.0 + vector1 += np.einsum('ia,a->i', cc.l1, cc.t1[p, :]) + vector1 += 2 * np.einsum('ilcd,lcd->i', cc.l2, cc.t2[p, :, :, :]) + vector1 -= np.einsum('ilcd,ldc->i', cc.l2, cc.t2[p, :, :, :]) + + vector2[p, :, :] += -2. * cc.l1 + vector2[:, p, :] += cc.l1 + vector2 += 2 * np.einsum('c,ijcb->ijb', cc.t1[p, :], cc.l2) + vector2 -= np.einsum('c,jicb->ijb', cc.t1[p, :], cc.l2) else: - vector1 += -cc.l1[:,p-nocc] - vector2 += -2*cc.l2[:,:,p-nocc,:] + cc.l2[:,:,:,p-nocc] + vector1 += -cc.l1[:, p - nocc] + vector2 += -2 * cc.l2[:, :, p - nocc, :] + cc.l2[:, :, :, p - nocc] - return amplitudes_to_vector_ip(vector1,vector2) + return amplitudes_to_vector_ip(vector1, vector2) -def greens_func_multiply(ham,vector,linear_part,args=None): - return np.array(ham(vector) + (linear_part)*vector) +def greens_func_multiply(ham, vector, linear_part, args=None): + return np.array(ham(vector) + (linear_part) * vector) def initial_ip_guess(cc): nocc, nvir = cc.t1.shape - vector1 = np.zeros((nocc),dtype=complex) - vector2 = np.zeros((nocc,nocc,nvir),dtype=complex) - return amplitudes_to_vector_ip(vector1,vector2) + vector1 = np.zeros((nocc), dtype=complex) + vector2 = np.zeros((nocc, nocc, nvir), dtype=complex) + return amplitudes_to_vector_ip(vector1, vector2) def initial_ea_guess(cc): nocc, nvir = cc.t1.shape - vector1 = np.zeros((nvir),dtype=complex) - vector2 = np.zeros((nocc,nvir,nvir),dtype=complex) - return amplitudes_to_vector_ea(vector1,vector2) + vector1 = np.zeros((nvir), dtype=complex) + vector2 = np.zeros((nocc, nvir, nvir), dtype=complex) + return amplitudes_to_vector_ea(vector1, vector2) def ip_shape(cc): nocc, nvir = cc.t1.shape - return nocc + nocc*nocc*nvir + return nocc + nocc * nocc * nvir def ea_shape(cc): nocc, nvir = cc.t1.shape - return nvir + nocc*nvir*nvir + return nvir + nocc * nvir * nvir class greens_function: def __init__(self, verbose=0): self.verbose = verbose - def td_ip_ao(self,cc,ps,times,mo_coeff,re_im="re",tol=1.e-5): + def td_ip_ao(self, cc, ps, times, mo_coeff, re_im="re", tol=1.e-5): """ E0: total CC gs energy ti: initial time @@ -150,35 +152,35 @@ def td_ip_ao(self,cc,ps,times,mo_coeff,re_im="re",tol=1.e-5): nmo = mo_coeff.shape[1] e_vector_mo = np.zeros([nmo, ip_shape(cc)], dtype=dtype) for i in range(nmo): - e_vector_mo[i,:] = greens_e_vector_ip_rhf(cc,i) - e_vector_ao = np.einsum("pi,ix->px", mo_coeff[ps,:], e_vector_mo) - + e_vector_mo[i, :] = greens_e_vector_ip_rhf(cc, i) + e_vector_ao = np.einsum("pi,ix->px", mo_coeff[ps, :], e_vector_mo) + b_vector_mo = np.zeros([ip_shape(cc), nmo], dtype=dtype) for i in range(nmo): - b_vector_mo[:,i] = greens_b_vector_ip_rhf(cc,i) - b_vector_ao = np.einsum("xi,ip->xp", b_vector_mo, mo_coeff.T[:,ps]) + b_vector_mo[:, i] = greens_b_vector_ip_rhf(cc, i) + b_vector_ao = np.einsum("xi,ip->xp", b_vector_mo, mo_coeff.T[:, ps]) # initialize loop variables - gf_ao = np.zeros((len(ps),len(ps),len(times)),dtype=dtype) - eomip=pyscf.cc.eom_rccsd.EOMIP(cc) - ti=times[0] - tf=times[-1] + gf_ao = np.zeros((len(ps), len(ps), len(times)), dtype=dtype) + eomip = pyscf.cc.eom_rccsd.EOMIP(cc) + ti = times[0] + tf = times[-1] - def matr_multiply(t,vector): + def matr_multiply(t, vector): # note: t is a dummy time argument, H is time-independent - return tfac*np.array(eomip.matvec(vector)) + return tfac * np.array(eomip.matvec(vector)) - for ip,p in enumerate(ps): - solp = scipy.integrate.solve_ivp(matr_multiply,(ti,tf), - b_vector_ao[:,p],t_eval=times, - rtol=tol,atol=tol) + for ip, p in enumerate(ps): + solp = scipy.integrate.solve_ivp(matr_multiply, (ti, tf), + b_vector_ao[:, p], t_eval=times, + rtol=tol, atol=tol) - for iq,q in enumerate(ps): - gf_ao[iq,ip,:] = np.dot(e_vector_ao[iq],solp.y) + for iq, q in enumerate(ps): + gf_ao[iq, ip, :] = np.dot(e_vector_ao[iq], solp.y) return gf_ao - def td_ea_ao(self,cc,ps,times,mo_coeff,re_im="re",tol=1.e-5): + def td_ea_ao(self, cc, ps, times, mo_coeff, re_im="re", tol=1.e-5): """ See td_ip. @@ -200,35 +202,35 @@ def td_ea_ao(self,cc,ps,times,mo_coeff,re_im="re",tol=1.e-5): nmo = mo_coeff.shape[1] e_vector_mo = np.zeros([nmo, ea_shape(cc)], dtype=dtype) for i in range(nmo): - e_vector_mo[i,:] = greens_e_vector_ea_rhf(cc,i) - e_vector_ao = np.einsum("pi,ix->px", mo_coeff[ps,:], e_vector_mo) - + e_vector_mo[i, :] = greens_e_vector_ea_rhf(cc, i) + e_vector_ao = np.einsum("pi,ix->px", mo_coeff[ps, :], e_vector_mo) + b_vector_mo = np.zeros([ea_shape(cc), nmo], dtype=dtype) for i in range(nmo): - b_vector_mo[:,i] = greens_b_vector_ea_rhf(cc,i) - b_vector_ao = np.einsum("xi,ip->xp", b_vector_mo, mo_coeff.T[:,ps]) + b_vector_mo[:, i] = greens_b_vector_ea_rhf(cc, i) + b_vector_ao = np.einsum("xi,ip->xp", b_vector_mo, mo_coeff.T[:, ps]) # initialize loop variables - gf_ao = np.zeros((len(ps),len(ps),len(times)),dtype=dtype) - eomea=pyscf.cc.eom_rccsd.EOMEA(cc) - ti=times[0] - tf=times[-1] + gf_ao = np.zeros((len(ps), len(ps), len(times)), dtype=dtype) + eomea = pyscf.cc.eom_rccsd.EOMEA(cc) + ti = times[0] + tf = times[-1] - def matr_multiply(t,vector): + def matr_multiply(t, vector): # note: t is a dummy time argument, H is time-independent - return tfac*np.array(eomea.matvec(vector)) + return tfac * np.array(eomea.matvec(vector)) - for ip,p in enumerate(ps): - solp = scipy.integrate.solve_ivp(matr_multiply,(ti,tf), - b_vector_ao[:,p],t_eval=times, - rtol=tol,atol=tol) + for ip, p in enumerate(ps): + solp = scipy.integrate.solve_ivp(matr_multiply, (ti, tf), + b_vector_ao[:, p], t_eval=times, + rtol=tol, atol=tol) - for iq,q in enumerate(ps): - gf_ao[iq,ip,:] = np.dot(e_vector_ao[iq],solp.y) + for iq, q in enumerate(ps): + gf_ao[iq, ip, :] = np.dot(e_vector_ao[iq], solp.y) return gf_ao - def td_ip(self,cc,ps,qs,times,re_im="re",tol=1.e-5): + def td_ip(self, cc, ps, qs, times, re_im="re", tol=1.e-5): """ E0: total CC gs energy ti: initial time @@ -252,34 +254,34 @@ def td_ip(self,cc,ps,qs,times,re_im="re",tol=1.e-5): tfac = 1j else: raise RuntimeError - - eomip=pyscf.cc.eom_rccsd.EOMIP(cc) - ti=times[0] - tf=times[-1] + + eomip = pyscf.cc.eom_rccsd.EOMIP(cc) + ti = times[0] + tf = times[-1] e_vector = list() for q in qs: - e_vector.append(greens_e_vector_ip_rhf(cc,q)) + e_vector.append(greens_e_vector_ip_rhf(cc, q)) - gfvals = np.zeros((len(ps),len(qs),len(times)),dtype=dtype) + gfvals = np.zeros((len(ps), len(qs), len(times)), dtype=dtype) - for ip,p in enumerate(ps): - b_vector = np.array(greens_b_vector_ip_rhf(cc,p), dtype=dtype) - - def matr_multiply(t,vector,args=None): + for ip, p in enumerate(ps): + b_vector = np.array(greens_b_vector_ip_rhf(cc, p), dtype=dtype) + + def matr_multiply(t, vector, args=None): # note: t is a dummy time argument, H is time-independent - res = tfac*np.array(eomip.matvec(vector)) + res = tfac * np.array(eomip.matvec(vector)) return res - - solp = scipy.integrate.solve_ivp(matr_multiply,(ti,tf), - b_vector,t_eval=times, rtol=tol, atol=tol) - for iq,q in enumerate(qs): - gfvals[iq,ip,:] = np.dot(e_vector[iq],solp.y) + solp = scipy.integrate.solve_ivp(matr_multiply, (ti, tf), + b_vector, t_eval=times, rtol=tol, atol=tol) + + for iq, q in enumerate(qs): + gfvals[iq, ip, :] = np.dot(e_vector[iq], solp.y) return gfvals - def td_ea(self,cc,ps,qs,times,re_im="re",tol=1.e-5): + def td_ea(self, cc, ps, qs, times, re_im="re", tol=1.e-5): """ See td_ip. @@ -297,145 +299,152 @@ def td_ea(self,cc,ps,qs,times,re_im="re",tol=1.e-5): else: raise RuntimeError - ti=times[0] - tf=times[-1] - eomea=pyscf.cc.eom_rccsd.EOMEA(cc) - + ti = times[0] + tf = times[-1] + eomea = pyscf.cc.eom_rccsd.EOMEA(cc) + e_vector = list() for p in ps: - e_vector.append(np.array(greens_e_vector_ea_rhf(cc,p), dtype=dtype)) - gfvals = np.zeros((len(ps),len(qs),len(times)),dtype=complex) - - for iq,q in enumerate(qs): - b_vector = np.array(greens_b_vector_ea_rhf(cc,q), dtype=dtype) + e_vector.append(np.array(greens_e_vector_ea_rhf(cc, p), dtype=dtype)) + gfvals = np.zeros((len(ps), len(qs), len(times)), dtype=complex) + + for iq, q in enumerate(qs): + b_vector = np.array(greens_b_vector_ea_rhf(cc, q), dtype=dtype) - def matr_multiply(t,vector,args=None): + def matr_multiply(t, vector, args=None): # t is a dummy time argument - res = tfac*np.array(eomea.matvec(vector)) + res = tfac * np.array(eomea.matvec(vector)) return res - solq = scipy.integrate.solve_ivp(matr_multiply,(ti,tf), - b_vector,t_eval=times, rtol=tol, atol=tol) - - for ip,p in enumerate(ps): - gfvals[ip,iq,:] = np.dot(e_vector[ip],solq.y) + solq = scipy.integrate.solve_ivp(matr_multiply, (ti, tf), + b_vector, t_eval=times, rtol=tol, atol=tol) + + for ip, p in enumerate(ps): + gfvals[ip, iq, :] = np.dot(e_vector[ip], solq.y) return gfvals - def solve_ip_ao(self,cc,ps,omega_list,mo_coeff,broadening): - eomip=pyscf.cc.eom_rccsd.EOMIP(cc) + def solve_ip_ao(self, cc, ps, omega_list, mo_coeff, broadening): + eomip = pyscf.cc.eom_rccsd.EOMIP(cc) # GKC: Why is this is the initial guess? x0 = initial_ip_guess(cc) - p0 = 0.0*x0 + 1.0 + p0 = 0.0 * x0 + 1.0 # set initial bra/ket nmo = mo_coeff.shape[1] e_vector_mo = np.zeros([nmo, ip_shape(cc)], dtype=np.complex128) for i in range(nmo): - e_vector_mo[i,:] = greens_e_vector_ip_rhf(cc,i) - e_vector_ao = np.einsum("pi,ix->px", mo_coeff[ps,:], e_vector_mo) + e_vector_mo[i, :] = greens_e_vector_ip_rhf(cc, i) + e_vector_ao = np.einsum("pi,ix->px", mo_coeff[ps, :], e_vector_mo) b_vector_mo = np.zeros([ip_shape(cc), nmo], dtype=np.complex128) for i in range(nmo): - b_vector_mo[:,i] = greens_b_vector_ip_rhf(cc,i) - b_vector_ao = np.einsum("xi,ip->xp", b_vector_mo, mo_coeff.T[:,ps]) + b_vector_mo[:, i] = greens_b_vector_ip_rhf(cc, i) + b_vector_ao = np.einsum("xi,ip->xp", b_vector_mo, mo_coeff.T[:, ps]) # initialize loop variables - gf_ao = np.zeros((len(ps),len(ps),len(omega_list)),dtype=np.complex128) - - for ip,p in enumerate(ps): + gf_ao = np.zeros((len(ps), len(ps), len(omega_list)), dtype=np.complex128) + + for ip, p in enumerate(ps): for iomega in range(len(omega_list)): curr_omega = omega_list[iomega] - def matr_multiply(vector,args=None): - return greens_func_multiply(eomip.matvec, vector, curr_omega-1j*broadening) - solver = gminres.gMinRes(matr_multiply,b_vector_ao[:,p],x0,p0) - #solver = gminres.exactInverse(matr_multiply,b_vector,x0) + + def matr_multiply(vector, args=None): + return greens_func_multiply(eomip.matvec, vector, curr_omega - 1j * broadening) + + solver = gminres.gMinRes(matr_multiply, b_vector_ao[:, p], x0, p0) + # solver = gminres.exactInverse(matr_multiply,b_vector,x0) sol = solver.get_solution().reshape(-1) - x0 = sol - for iq,q in enumerate(ps): - gf_ao[ip,iq,iomega] = -np.dot(e_vector_ao[iq,:],sol) + x0 = sol + for iq, q in enumerate(ps): + gf_ao[ip, iq, iomega] = -np.dot(e_vector_ao[iq, :], sol) return gf_ao - def solve_ea_ao(self,cc,ps,omega_list,mo_coeff,broadening): - eomea=pyscf.cc.eom_rccsd.EOMEA(cc) + def solve_ea_ao(self, cc, ps, omega_list, mo_coeff, broadening): + eomea = pyscf.cc.eom_rccsd.EOMEA(cc) # GKC: Why is this is the initial guess? x0 = initial_ea_guess(cc) - p0 = 0.0*x0 + 1.0 + p0 = 0.0 * x0 + 1.0 # set initial bra/ket nmo = mo_coeff.shape[1] e_vector_mo = np.zeros([nmo, ea_shape(cc)], dtype=np.complex128) for i in range(nmo): - e_vector_mo[i,:] = greens_e_vector_ea_rhf(cc,i) - e_vector_ao = np.einsum("pi,ix->px", mo_coeff[ps,:], e_vector_mo) + e_vector_mo[i, :] = greens_e_vector_ea_rhf(cc, i) + e_vector_ao = np.einsum("pi,ix->px", mo_coeff[ps, :], e_vector_mo) b_vector_mo = np.zeros([ea_shape(cc), nmo], dtype=np.complex128) for i in range(nmo): - b_vector_mo[:,i] = greens_b_vector_ea_rhf(cc,i) - b_vector_ao = np.einsum("xi,ip->xp", b_vector_mo, mo_coeff.T[:,ps]) + b_vector_mo[:, i] = greens_b_vector_ea_rhf(cc, i) + b_vector_ao = np.einsum("xi,ip->xp", b_vector_mo, mo_coeff.T[:, ps]) # initialize loop variables - gf_ao = np.zeros((len(ps),len(ps),len(omega_list)),dtype=np.complex128) - - for iq,q in enumerate(ps): + gf_ao = np.zeros((len(ps), len(ps), len(omega_list)), dtype=np.complex128) + + for iq, q in enumerate(ps): for iomega in range(len(omega_list)): curr_omega = omega_list[iomega] - def matr_multiply(vector,args=None): - return greens_func_multiply(eomea.matvec, vector, -curr_omega-1j*broadening) - solver = gminres.gMinRes(matr_multiply,b_vector_ao[:,q],x0,p0) - #solver = gminres.exactInverse(matr_multiply,b_vector,x0) + + def matr_multiply(vector, args=None): + return greens_func_multiply(eomea.matvec, vector, -curr_omega - 1j * broadening) + + solver = gminres.gMinRes(matr_multiply, b_vector_ao[:, q], x0, p0) + # solver = gminres.exactInverse(matr_multiply,b_vector,x0) sol = solver.get_solution().reshape(-1) - x0 = sol - for ip,p in enumerate(ps): - gf_ao[ip,iq,iomega] = np.dot(e_vector_ao[ip],sol) + x0 = sol + for ip, p in enumerate(ps): + gf_ao[ip, iq, iomega] = np.dot(e_vector_ao[ip], sol) return gf_ao - - def solve_ip(self,cc,ps,qs,omega_list,broadening): - eomip=pyscf.cc.eom_rccsd.EOMIP(cc) + + def solve_ip(self, cc, ps, qs, omega_list, broadening): + eomip = pyscf.cc.eom_rccsd.EOMIP(cc) x0 = initial_ip_guess(cc) - p0 = 0.0*x0 + 1.0 - e_vector = list() + p0 = 0.0 * x0 + 1.0 + e_vector = list() for q in qs: - e_vector.append(greens_e_vector_ip_rhf(cc,q)) - gfvals = np.zeros((len(ps),len(qs),len(omega_list)),dtype=complex) - for ip,p in enumerate(ps): - b_vector = greens_b_vector_ip_rhf(cc,p) + e_vector.append(greens_e_vector_ip_rhf(cc, q)) + gfvals = np.zeros((len(ps), len(qs), len(omega_list)), dtype=complex) + for ip, p in enumerate(ps): + b_vector = greens_b_vector_ip_rhf(cc, p) for iomega in range(len(omega_list)): curr_omega = omega_list[iomega] - def matr_multiply(vector,args=None): - return greens_func_multiply(eomip.matvec, vector, curr_omega-1j*broadening) - solver = gminres.gMinRes(matr_multiply,b_vector,x0,p0) - #solver = gminres.exactInverse(matr_multiply,b_vector,x0) + + def matr_multiply(vector, args=None): + return greens_func_multiply(eomip.matvec, vector, curr_omega - 1j * broadening) + + solver = gminres.gMinRes(matr_multiply, b_vector, x0, p0) + # solver = gminres.exactInverse(matr_multiply,b_vector,x0) sol = solver.get_solution().reshape(-1) - x0 = sol - for iq,q in enumerate(qs): - gfvals[ip,iq,iomega] = -np.dot(e_vector[iq],sol) + x0 = sol + for iq, q in enumerate(qs): + gfvals[ip, iq, iomega] = -np.dot(e_vector[iq], sol) if len(ps) == 1 and len(qs) == 1: - return gfvals[0,0,:] + return gfvals[0, 0, :] else: return gfvals - def solve_ea(self,cc,ps,qs,omega_list,broadening): - eomea=pyscf.cc.eom_rccsd.EOMEA(cc) + def solve_ea(self, cc, ps, qs, omega_list, broadening): + eomea = pyscf.cc.eom_rccsd.EOMEA(cc) x0 = initial_ea_guess(cc) - p0 = 0.0*x0 + 1.0 - e_vector = list() + p0 = 0.0 * x0 + 1.0 + e_vector = list() for p in ps: - e_vector.append(greens_e_vector_ea_rhf(cc,p)) - gfvals = np.zeros((len(ps),len(qs),len(omega_list)),dtype=complex) - for iq,q in enumerate(qs): - b_vector = greens_b_vector_ea_rhf(cc,q) + e_vector.append(greens_e_vector_ea_rhf(cc, p)) + gfvals = np.zeros((len(ps), len(qs), len(omega_list)), dtype=complex) + for iq, q in enumerate(qs): + b_vector = greens_b_vector_ea_rhf(cc, q) for iomega in range(len(omega_list)): curr_omega = omega_list[iomega] - def matr_multiply(vector,args=None): - return greens_func_multiply(eomea.matvec, vector, -curr_omega-1j*broadening) - solver = gminres.gMinRes(matr_multiply,b_vector,x0,p0) - #solver = gminres.exactInverse(matr_multiply,b_vector,x0) + + def matr_multiply(vector, args=None): + return greens_func_multiply(eomea.matvec, vector, -curr_omega - 1j * broadening) + + solver = gminres.gMinRes(matr_multiply, b_vector, x0, p0) + # solver = gminres.exactInverse(matr_multiply,b_vector,x0) sol = solver.get_solution().reshape(-1) x0 = sol - for ip,p in enumerate(ps): - gfvals[ip,iq,iomega] = np.dot(e_vector[ip],sol) + for ip, p in enumerate(ps): + gfvals[ip, iq, iomega] = np.dot(e_vector[ip], sol) if len(ps) == 1 and len(qs) == 1: - return gfvals[0,0,:] + return gfvals[0, 0, :] else: return gfvals - def solve_gf(self,cc,p,q,omega_list,broadening): - return self.solve_ip(cc,p,q,omega_list,broadening), self.solve_ea(cc,p,q,omega_list,broadening) - + def solve_gf(self, cc, p, q, omega_list, broadening): + return self.solve_ip(cc, p, q, omega_list, broadening), self.solve_ea(cc, p, q, omega_list, broadening) From 7e600bfaad1a2540ad481fc53a0cbaa276d76a9d Mon Sep 17 00:00:00 2001 From: Unknown Date: Tue, 24 Jul 2018 17:56:59 -0700 Subject: [PATCH 3/7] Missing comma typos fix --- greens_function.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/greens_function.py b/greens_function.py index d37cce5..e6440f9 100644 --- a/greens_function.py +++ b/greens_function.py @@ -15,7 +15,7 @@ def greens_b_vector_ea_rhf(cc, p): nocc, nvir = cc.t1.shape ds_type = cc.t1.dtype - vector1 = np.zeros((nvir), dtype=ds_type) + vector1 = np.zeros((nvir,), dtype=ds_type) vector2 = np.zeros((nocc, nvir, nvir), dtype=ds_type) if p < nocc: # Changed both to minus @@ -29,7 +29,7 @@ def greens_b_vector_ea_rhf(cc, p): def greens_e_vector_ea_rhf(cc, p): nocc, nvir = cc.t1.shape ds_type = cc.t1.dtype - vector1 = np.zeros((nvir), dtype=ds_type) + vector1 = np.zeros((nvir,), dtype=ds_type) vector2 = np.zeros((nocc, nvir, nvir), dtype=ds_type) if p < nocc: # Changed both to plus @@ -56,7 +56,7 @@ def greens_e_vector_ea_rhf(cc, p): def greens_b_vector_ip_rhf(cc, p): nocc, nvir = cc.t1.shape - vector1 = np.zeros((nocc), dtype=complex) + vector1 = np.zeros((nocc,), dtype=complex) vector2 = np.zeros((nocc, nocc, nvir), dtype=complex) if p < nocc: vector1[p] = 1.0 @@ -68,7 +68,7 @@ def greens_b_vector_ip_rhf(cc, p): def greens_e_vector_ip_rhf(cc, p): nocc, nvir = cc.t1.shape - vector1 = np.zeros((nocc), dtype=complex) + vector1 = np.zeros((nocc,), dtype=complex) vector2 = np.zeros((nocc, nocc, nvir), dtype=complex) if p < nocc: vector1[p] = -1.0 @@ -88,19 +88,19 @@ def greens_e_vector_ip_rhf(cc, p): def greens_func_multiply(ham, vector, linear_part, args=None): - return np.array(ham(vector) + (linear_part) * vector) + return np.array(ham(vector) + linear_part * vector) def initial_ip_guess(cc): nocc, nvir = cc.t1.shape - vector1 = np.zeros((nocc), dtype=complex) + vector1 = np.zeros((nocc,), dtype=complex) vector2 = np.zeros((nocc, nocc, nvir), dtype=complex) return amplitudes_to_vector_ip(vector1, vector2) def initial_ea_guess(cc): nocc, nvir = cc.t1.shape - vector1 = np.zeros((nvir), dtype=complex) + vector1 = np.zeros((nvir,), dtype=complex) vector2 = np.zeros((nocc, nvir, nvir), dtype=complex) return amplitudes_to_vector_ea(vector1, vector2) From b3032064c035eb32538c0c8ddb8cb63554c1f300 Mon Sep 17 00:00:00 2001 From: Unknown Date: Tue, 24 Jul 2018 18:43:43 -0700 Subject: [PATCH 4/7] Removed code repetition and added an interface to scipy linear solver --- greens_function.py | 42 +++++++++++++++++++++++++++++------------- 1 file changed, 29 insertions(+), 13 deletions(-) diff --git a/greens_function.py b/greens_function.py index e6440f9..6cec254 100644 --- a/greens_function.py +++ b/greens_function.py @@ -1,6 +1,5 @@ import numpy as np import scipy -import gminres import pyscf import pyscf.cc @@ -349,9 +348,7 @@ def solve_ip_ao(self, cc, ps, omega_list, mo_coeff, broadening): def matr_multiply(vector, args=None): return greens_func_multiply(eomip.matvec, vector, curr_omega - 1j * broadening) - solver = gminres.gMinRes(matr_multiply, b_vector_ao[:, p], x0, p0) - # solver = gminres.exactInverse(matr_multiply,b_vector,x0) - sol = solver.get_solution().reshape(-1) + sol = self.__solve_linear_problem__(matr_multiply, b_vector_ao[:, p], x0, p0) x0 = sol for iq, q in enumerate(ps): gf_ao[ip, iq, iomega] = -np.dot(e_vector_ao[iq, :], sol) @@ -383,9 +380,7 @@ def solve_ea_ao(self, cc, ps, omega_list, mo_coeff, broadening): def matr_multiply(vector, args=None): return greens_func_multiply(eomea.matvec, vector, -curr_omega - 1j * broadening) - solver = gminres.gMinRes(matr_multiply, b_vector_ao[:, q], x0, p0) - # solver = gminres.exactInverse(matr_multiply,b_vector,x0) - sol = solver.get_solution().reshape(-1) + sol = self.__solve_linear_problem__(matr_multiply, b_vector_ao[:, q], x0, p0) x0 = sol for ip, p in enumerate(ps): gf_ao[ip, iq, iomega] = np.dot(e_vector_ao[ip], sol) @@ -408,9 +403,7 @@ def solve_ip(self, cc, ps, qs, omega_list, broadening): def matr_multiply(vector, args=None): return greens_func_multiply(eomip.matvec, vector, curr_omega - 1j * broadening) - solver = gminres.gMinRes(matr_multiply, b_vector, x0, p0) - # solver = gminres.exactInverse(matr_multiply,b_vector,x0) - sol = solver.get_solution().reshape(-1) + sol = self.__solve_linear_problem__(matr_multiply, b_vector, x0, p0) x0 = sol for iq, q in enumerate(qs): gfvals[ip, iq, iomega] = -np.dot(e_vector[iq], sol) @@ -435,9 +428,7 @@ def solve_ea(self, cc, ps, qs, omega_list, broadening): def matr_multiply(vector, args=None): return greens_func_multiply(eomea.matvec, vector, -curr_omega - 1j * broadening) - solver = gminres.gMinRes(matr_multiply, b_vector, x0, p0) - # solver = gminres.exactInverse(matr_multiply,b_vector,x0) - sol = solver.get_solution().reshape(-1) + sol = self.__solve_linear_problem__(matr_multiply, b_vector, x0, p0) x0 = sol for ip, p in enumerate(ps): gfvals[ip, iq, iomega] = np.dot(e_vector[ip], sol) @@ -448,3 +439,28 @@ def matr_multiply(vector, args=None): def solve_gf(self, cc, p, q, omega_list, broadening): return self.solve_ip(cc, p, q, omega_list, broadening), self.solve_ea(cc, p, q, omega_list, broadening) + + def __solve_linear_problem__(self, matr_multiply, b_vector, x0, p0): + self.l(" Solving linear problem ...") + try: + import gminres + driver = gminres.gMinRes(matr_multiply, b_vector, x0, p0) + return driver.get_solution().reshape(-1) + + except ImportError: + self.l(" gminres import failed, trying scipy ...") + # TODO: p0 is not used here, whatever it is + import scipy.sparse.linalg as spla + size = len(b_vector) + Ax = spla.LinearOperator((size, size), matr_multiply) + result, info = spla.gmres(Ax, b_vector, x0=x0, callback=self.__spla_callback__) + if info != 0: + raise RuntimeError("Error solving linear problem: info={:d}".format(info)) + return result + + def l(self, x): + if self.verbose > 0: + print(x) + + def __spla_callback__(self, r): + self.l(" res = {:.3e}".format(r)) \ No newline at end of file From 2efa2987fe933eb60182a74f0d74e4d02578fcc3 Mon Sep 17 00:00:00 2001 From: Unknown Date: Wed, 25 Jul 2018 17:02:18 -0700 Subject: [PATCH 5/7] Cached intermediates for a speedup --- greens_function.py | 28 ++++++++++++++++++---------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/greens_function.py b/greens_function.py index 6cec254..7e5e3f1 100644 --- a/greens_function.py +++ b/greens_function.py @@ -86,8 +86,8 @@ def greens_e_vector_ip_rhf(cc, p): return amplitudes_to_vector_ip(vector1, vector2) -def greens_func_multiply(ham, vector, linear_part, args=None): - return np.array(ham(vector) + linear_part * vector) +def greens_func_multiply(ham, vector, linear_part, **kwargs): + return np.array(ham(vector, **kwargs) + linear_part * vector) def initial_ip_guess(cc): @@ -162,12 +162,13 @@ def td_ip_ao(self, cc, ps, times, mo_coeff, re_im="re", tol=1.e-5): # initialize loop variables gf_ao = np.zeros((len(ps), len(ps), len(times)), dtype=dtype) eomip = pyscf.cc.eom_rccsd.EOMIP(cc) + eomip_imds = eomip.make_imds() ti = times[0] tf = times[-1] def matr_multiply(t, vector): # note: t is a dummy time argument, H is time-independent - return tfac * np.array(eomip.matvec(vector)) + return tfac * np.array(eomip.matvec(vector, imds=eomip_imds)) for ip, p in enumerate(ps): solp = scipy.integrate.solve_ivp(matr_multiply, (ti, tf), @@ -212,12 +213,13 @@ def td_ea_ao(self, cc, ps, times, mo_coeff, re_im="re", tol=1.e-5): # initialize loop variables gf_ao = np.zeros((len(ps), len(ps), len(times)), dtype=dtype) eomea = pyscf.cc.eom_rccsd.EOMEA(cc) + eomea_imds = eomea.make_imds() ti = times[0] tf = times[-1] def matr_multiply(t, vector): # note: t is a dummy time argument, H is time-independent - return tfac * np.array(eomea.matvec(vector)) + return tfac * np.array(eomea.matvec(vector, imds=eomea_imds)) for ip, p in enumerate(ps): solp = scipy.integrate.solve_ivp(matr_multiply, (ti, tf), @@ -255,6 +257,7 @@ def td_ip(self, cc, ps, qs, times, re_im="re", tol=1.e-5): raise RuntimeError eomip = pyscf.cc.eom_rccsd.EOMIP(cc) + eomip_imds = eomip.make_imds() ti = times[0] tf = times[-1] @@ -269,7 +272,7 @@ def td_ip(self, cc, ps, qs, times, re_im="re", tol=1.e-5): def matr_multiply(t, vector, args=None): # note: t is a dummy time argument, H is time-independent - res = tfac * np.array(eomip.matvec(vector)) + res = tfac * np.array(eomip.matvec(vector, imds=eomip_imds)) return res solp = scipy.integrate.solve_ivp(matr_multiply, (ti, tf), @@ -301,6 +304,7 @@ def td_ea(self, cc, ps, qs, times, re_im="re", tol=1.e-5): ti = times[0] tf = times[-1] eomea = pyscf.cc.eom_rccsd.EOMEA(cc) + eomea_imds = eomea.make_imds() e_vector = list() for p in ps: @@ -312,7 +316,7 @@ def td_ea(self, cc, ps, qs, times, re_im="re", tol=1.e-5): def matr_multiply(t, vector, args=None): # t is a dummy time argument - res = tfac * np.array(eomea.matvec(vector)) + res = tfac * np.array(eomea.matvec(vector, imds=eomea_imds)) return res solq = scipy.integrate.solve_ivp(matr_multiply, (ti, tf), @@ -324,6 +328,7 @@ def matr_multiply(t, vector, args=None): def solve_ip_ao(self, cc, ps, omega_list, mo_coeff, broadening): eomip = pyscf.cc.eom_rccsd.EOMIP(cc) + eomip_imds = eomip.make_imds() # GKC: Why is this is the initial guess? x0 = initial_ip_guess(cc) p0 = 0.0 * x0 + 1.0 @@ -346,7 +351,7 @@ def solve_ip_ao(self, cc, ps, omega_list, mo_coeff, broadening): curr_omega = omega_list[iomega] def matr_multiply(vector, args=None): - return greens_func_multiply(eomip.matvec, vector, curr_omega - 1j * broadening) + return greens_func_multiply(eomip.matvec, vector, curr_omega - 1j * broadening, imds=eomip_imds) sol = self.__solve_linear_problem__(matr_multiply, b_vector_ao[:, p], x0, p0) x0 = sol @@ -356,6 +361,7 @@ def matr_multiply(vector, args=None): def solve_ea_ao(self, cc, ps, omega_list, mo_coeff, broadening): eomea = pyscf.cc.eom_rccsd.EOMEA(cc) + eomea_imds = eomea.make_imds() # GKC: Why is this is the initial guess? x0 = initial_ea_guess(cc) p0 = 0.0 * x0 + 1.0 @@ -378,7 +384,7 @@ def solve_ea_ao(self, cc, ps, omega_list, mo_coeff, broadening): curr_omega = omega_list[iomega] def matr_multiply(vector, args=None): - return greens_func_multiply(eomea.matvec, vector, -curr_omega - 1j * broadening) + return greens_func_multiply(eomea.matvec, vector, -curr_omega - 1j * broadening, imds=eomea_imds) sol = self.__solve_linear_problem__(matr_multiply, b_vector_ao[:, q], x0, p0) x0 = sol @@ -389,6 +395,7 @@ def matr_multiply(vector, args=None): def solve_ip(self, cc, ps, qs, omega_list, broadening): eomip = pyscf.cc.eom_rccsd.EOMIP(cc) + eomip_imds = eomip.make_imds() x0 = initial_ip_guess(cc) p0 = 0.0 * x0 + 1.0 e_vector = list() @@ -401,7 +408,7 @@ def solve_ip(self, cc, ps, qs, omega_list, broadening): curr_omega = omega_list[iomega] def matr_multiply(vector, args=None): - return greens_func_multiply(eomip.matvec, vector, curr_omega - 1j * broadening) + return greens_func_multiply(eomip.matvec, vector, curr_omega - 1j * broadening, imds=eomip_imds) sol = self.__solve_linear_problem__(matr_multiply, b_vector, x0, p0) x0 = sol @@ -414,6 +421,7 @@ def matr_multiply(vector, args=None): def solve_ea(self, cc, ps, qs, omega_list, broadening): eomea = pyscf.cc.eom_rccsd.EOMEA(cc) + eomea_imds = eomea.make_imds() x0 = initial_ea_guess(cc) p0 = 0.0 * x0 + 1.0 e_vector = list() @@ -426,7 +434,7 @@ def solve_ea(self, cc, ps, qs, omega_list, broadening): curr_omega = omega_list[iomega] def matr_multiply(vector, args=None): - return greens_func_multiply(eomea.matvec, vector, -curr_omega - 1j * broadening) + return greens_func_multiply(eomea.matvec, vector, -curr_omega - 1j * broadening, imds=eomea_imds) sol = self.__solve_linear_problem__(matr_multiply, b_vector, x0, p0) x0 = sol From ebcc7ff9d6dfd899bd019417c8cf935f4eed776f Mon Sep 17 00:00:00 2001 From: Unknown Date: Thu, 26 Jul 2018 16:01:40 -0700 Subject: [PATCH 6/7] Split logic into numpy part and pyscf part for vector functions --- greens_function.py | 110 +++++++++++++++++++++++++-------------------- 1 file changed, 61 insertions(+), 49 deletions(-) diff --git a/greens_function.py b/greens_function.py index 7e5e3f1..d2476a2 100644 --- a/greens_function.py +++ b/greens_function.py @@ -11,41 +11,45 @@ ################### -def greens_b_vector_ea_rhf(cc, p): - nocc, nvir = cc.t1.shape - ds_type = cc.t1.dtype - vector1 = np.zeros((nvir,), dtype=ds_type) - vector2 = np.zeros((nocc, nvir, nvir), dtype=ds_type) +def greens_b_amplitudes_ea_rhf(t1, t2, p): + nocc, nvir = t1.shape + ds_type = t1.dtype if p < nocc: - # Changed both to minus - vector1 += -cc.t1[p, :] - vector2 += -cc.t2[p, :, :, :] + return -t1[p, :], -t2[p, :, :, :] else: + vector1 = np.zeros((nvir,), dtype=ds_type) + vector2 = np.zeros((nocc, nvir, nvir), dtype=ds_type) vector1[p - nocc] = 1.0 - return amplitudes_to_vector_ea(vector1, vector2) + return vector1, vector2 -def greens_e_vector_ea_rhf(cc, p): - nocc, nvir = cc.t1.shape - ds_type = cc.t1.dtype - vector1 = np.zeros((nvir,), dtype=ds_type) - vector2 = np.zeros((nocc, nvir, nvir), dtype=ds_type) +def greens_b_vector_ea_rhf(cc, p): + return amplitudes_to_vector_ea(*greens_b_amplitudes_ea_rhf(cc.t1, cc.t2, p)) + + +def greens_e_amplitudes_ea_rhf(t1, t2, l1, l2, p): + nocc, nvir = t1.shape + ds_type = t1.dtype if p < nocc: - # Changed both to plus - vector1 += cc.l1[p, :] - vector2 += (2 * cc.l2[p, :, :, :] - cc.l2[:, p, :, :]) - pass + return l1[p, :], 2 * l2[p, :, :, :] - l2[:, p, :, :] else: + vector1 = np.zeros((nvir,), dtype=ds_type) + vector2 = np.zeros((nocc, nvir, nvir), dtype=ds_type) vector1[p - nocc] = -1.0 - vector1 += np.einsum('ia,i->a', cc.l1, cc.t1[:, p - nocc]) - vector1 += 2 * np.einsum('klca,klc->a', cc.l2, cc.t2[:, :, :, p - nocc]) - vector1 -= np.einsum('klca,lkc->a', cc.l2, cc.t2[:, :, :, p - nocc]) - - vector2[:, p - nocc, :] += -2. * cc.l1 - vector2[:, :, p - nocc] += cc.l1 - vector2 += 2 * np.einsum('k,jkba->jab', cc.t1[:, p - nocc], cc.l2) - vector2 -= np.einsum('k,jkab->jab', cc.t1[:, p - nocc], cc.l2) - return amplitudes_to_vector_ea(vector1, vector2) + vector1 += np.einsum('ia,i->a', l1, t1[:, p - nocc]) + vector1 += 2 * np.einsum('klca,klc->a', l2, t2[:, :, :, p - nocc]) + vector1 -= np.einsum('klca,lkc->a', l2, t2[:, :, :, p - nocc]) + + vector2[:, p - nocc, :] += -2. * l1 + vector2[:, :, p - nocc] += l1 + vector2 += 2 * np.einsum('k,jkba->jab', t1[:, p - nocc], l2) + vector2 -= np.einsum('k,jkab->jab', t1[:, p - nocc], l2) + + return vector1, vector2 + + +def greens_e_vector_ea_rhf(cc, p): + return amplitudes_to_vector_ea(*greens_e_amplitudes_ea_rhf(cc.t1, cc.t2, cc.l1, cc.l2, p)) ################### @@ -53,37 +57,45 @@ def greens_e_vector_ea_rhf(cc, p): ################### -def greens_b_vector_ip_rhf(cc, p): - nocc, nvir = cc.t1.shape - vector1 = np.zeros((nocc,), dtype=complex) - vector2 = np.zeros((nocc, nocc, nvir), dtype=complex) +def greens_b_amplitudes_ip_rhf(t1, t2, p): + nocc, nvir = t1.shape + ds_type = t1.dtype if p < nocc: + vector1 = np.zeros((nocc,), dtype=ds_type) + vector2 = np.zeros((nocc, nocc, nvir), dtype=ds_type) vector1[p] = 1.0 + return vector1, vector2 else: - vector1 += cc.t1[:, p - nocc] - vector2 += cc.t2[:, :, :, p - nocc] - return amplitudes_to_vector_ip(vector1, vector2) + return t1[:, p - nocc], t2[:, :, :, p - nocc] -def greens_e_vector_ip_rhf(cc, p): - nocc, nvir = cc.t1.shape - vector1 = np.zeros((nocc,), dtype=complex) - vector2 = np.zeros((nocc, nocc, nvir), dtype=complex) +def greens_b_vector_ip_rhf(cc, p): + return amplitudes_to_vector_ip(*greens_b_amplitudes_ip_rhf(cc.t1, cc.t2, p)) + + +def greens_e_amplitudes_ip_rhf(t1, t2, l1, l2, p): + nocc, nvir = t1.shape + ds_type = t1.dtype if p < nocc: + vector1 = np.zeros((nocc,), dtype=ds_type) + vector2 = np.zeros((nocc, nocc, nvir), dtype=ds_type) vector1[p] = -1.0 - vector1 += np.einsum('ia,a->i', cc.l1, cc.t1[p, :]) - vector1 += 2 * np.einsum('ilcd,lcd->i', cc.l2, cc.t2[p, :, :, :]) - vector1 -= np.einsum('ilcd,ldc->i', cc.l2, cc.t2[p, :, :, :]) - - vector2[p, :, :] += -2. * cc.l1 - vector2[:, p, :] += cc.l1 - vector2 += 2 * np.einsum('c,ijcb->ijb', cc.t1[p, :], cc.l2) - vector2 -= np.einsum('c,jicb->ijb', cc.t1[p, :], cc.l2) + vector1 += np.einsum('ia,a->i', l1, t1[p, :]) + vector1 += 2 * np.einsum('ilcd,lcd->i', l2, t2[p, :, :, :]) + vector1 -= np.einsum('ilcd,ldc->i', l2, t2[p, :, :, :]) + + vector2[p, :, :] += -2. * l1 + vector2[:, p, :] += l1 + vector2 += 2 * np.einsum('c,ijcb->ijb', t1[p, :], l2) + vector2 -= np.einsum('c,jicb->ijb', t1[p, :], l2) + + return vector1, vector2 else: - vector1 += -cc.l1[:, p - nocc] - vector2 += -2 * cc.l2[:, :, p - nocc, :] + cc.l2[:, :, :, p - nocc] + return -l1[:, p - nocc], -2 * l2[:, :, p - nocc, :] + l2[:, :, :, p - nocc] - return amplitudes_to_vector_ip(vector1, vector2) + +def greens_e_vector_ip_rhf(cc, p): + return amplitudes_to_vector_ip(*greens_e_amplitudes_ip_rhf(cc.t1, cc.t2, cc.l1, cc.l2, p)) def greens_func_multiply(ham, vector, linear_part, **kwargs): From 60727fb0ad9537ca86f948989c873dfe3f03524a Mon Sep 17 00:00:00 2001 From: Unknown Date: Thu, 26 Jul 2018 16:23:52 -0700 Subject: [PATCH 7/7] Further decomposed functions into singles and doubles parts --- greens_function.py | 130 ++++++++++++++++++++++++++++++--------------- 1 file changed, 86 insertions(+), 44 deletions(-) diff --git a/greens_function.py b/greens_function.py index d2476a2..9cc3b15 100644 --- a/greens_function.py +++ b/greens_function.py @@ -11,45 +11,66 @@ ################### -def greens_b_amplitudes_ea_rhf(t1, t2, p): +def greens_b_singles_ea_rhf(t1, p): nocc, nvir = t1.shape ds_type = t1.dtype if p < nocc: - return -t1[p, :], -t2[p, :, :, :] + return -t1[p, :] else: - vector1 = np.zeros((nvir,), dtype=ds_type) - vector2 = np.zeros((nocc, nvir, nvir), dtype=ds_type) - vector1[p - nocc] = 1.0 - return vector1, vector2 + result = np.zeros((nvir,), dtype=ds_type) + result[p - nocc] = 1.0 + return result + + +def greens_b_doubles_ea_rhf(t2, p): + nocc, _, nvir, _ = t2.shape + ds_type = t2.dtype + if p < nocc: + return -t2[p, :, :, :] + else: + return np.zeros((nocc, nvir, nvir), dtype=ds_type) def greens_b_vector_ea_rhf(cc, p): - return amplitudes_to_vector_ea(*greens_b_amplitudes_ea_rhf(cc.t1, cc.t2, p)) + return amplitudes_to_vector_ea( + greens_b_singles_ea_rhf(cc.t1, p), + greens_b_doubles_ea_rhf(cc.t2, p), + ) -def greens_e_amplitudes_ea_rhf(t1, t2, l1, l2, p): +def greens_e_singles_ea_rhf(t1, t2, l1, l2, p): nocc, nvir = t1.shape ds_type = t1.dtype if p < nocc: - return l1[p, :], 2 * l2[p, :, :, :] - l2[:, p, :, :] + return l1[p, :] else: - vector1 = np.zeros((nvir,), dtype=ds_type) - vector2 = np.zeros((nocc, nvir, nvir), dtype=ds_type) - vector1[p - nocc] = -1.0 - vector1 += np.einsum('ia,i->a', l1, t1[:, p - nocc]) - vector1 += 2 * np.einsum('klca,klc->a', l2, t2[:, :, :, p - nocc]) - vector1 -= np.einsum('klca,lkc->a', l2, t2[:, :, :, p - nocc]) + result = np.zeros((nvir,), dtype=ds_type) + result[p - nocc] = -1.0 + result += np.einsum('ia,i->a', l1, t1[:, p - nocc]) + result += 2 * np.einsum('klca,klc->a', l2, t2[:, :, :, p - nocc]) + result -= np.einsum('klca,lkc->a', l2, t2[:, :, :, p - nocc]) + return result - vector2[:, p - nocc, :] += -2. * l1 - vector2[:, :, p - nocc] += l1 - vector2 += 2 * np.einsum('k,jkba->jab', t1[:, p - nocc], l2) - vector2 -= np.einsum('k,jkab->jab', t1[:, p - nocc], l2) - return vector1, vector2 +def greens_e_doubles_ea_rhf(t1, l1, l2, p): + nocc, nvir = t1.shape + ds_type = t1.dtype + if p < nocc: + return 2 * l2[p, :, :, :] - l2[:, p, :, :] + else: + result = np.zeros((nocc, nvir, nvir), dtype=ds_type) + result[:, p - nocc, :] += -2. * l1 + result[:, :, p - nocc] += l1 + result += 2 * np.einsum('k,jkba->jab', t1[:, p - nocc], l2) + result -= np.einsum('k,jkab->jab', t1[:, p - nocc], l2) + return result def greens_e_vector_ea_rhf(cc, p): - return amplitudes_to_vector_ea(*greens_e_amplitudes_ea_rhf(cc.t1, cc.t2, cc.l1, cc.l2, p)) + return amplitudes_to_vector_ea( + greens_e_singles_ea_rhf(cc.t1, cc.t2, cc.l1, cc.l2, p), + greens_e_doubles_ea_rhf(cc.t1, cc.l1, cc.l2, p), + ) ################### @@ -57,45 +78,66 @@ def greens_e_vector_ea_rhf(cc, p): ################### -def greens_b_amplitudes_ip_rhf(t1, t2, p): +def greens_b_singles_ip_rhf(t1, p): nocc, nvir = t1.shape ds_type = t1.dtype if p < nocc: - vector1 = np.zeros((nocc,), dtype=ds_type) - vector2 = np.zeros((nocc, nocc, nvir), dtype=ds_type) - vector1[p] = 1.0 - return vector1, vector2 + result = np.zeros((nocc,), dtype=ds_type) + result[p] = 1.0 + return result + else: + return t1[:, p - nocc] + + +def greens_b_doubles_ip_rhf(t2, p): + nocc, _, nvir, _ = t2.shape + ds_type = t2.dtype + if p < nocc: + return np.zeros((nocc, nocc, nvir), dtype=ds_type) else: - return t1[:, p - nocc], t2[:, :, :, p - nocc] + return t2[:, :, :, p - nocc] def greens_b_vector_ip_rhf(cc, p): - return amplitudes_to_vector_ip(*greens_b_amplitudes_ip_rhf(cc.t1, cc.t2, p)) + return amplitudes_to_vector_ip( + greens_b_singles_ip_rhf(cc.t1, p), + greens_b_doubles_ip_rhf(cc.t2, p), + ) + + +def greens_e_singles_ip_rhf(t1, t2, l1, l2, p): + nocc, nvir = t1.shape + ds_type = t1.dtype + if p < nocc: + result = np.zeros((nocc,), dtype=ds_type) + result[p] = -1.0 + result += np.einsum('ia,a->i', l1, t1[p, :]) + result += 2 * np.einsum('ilcd,lcd->i', l2, t2[p, :, :, :]) + result -= np.einsum('ilcd,ldc->i', l2, t2[p, :, :, :]) + return result + else: + return -l1[:, p - nocc] -def greens_e_amplitudes_ip_rhf(t1, t2, l1, l2, p): +def greens_e_doubles_ip_rhf(t1, l1, l2, p): nocc, nvir = t1.shape ds_type = t1.dtype if p < nocc: - vector1 = np.zeros((nocc,), dtype=ds_type) - vector2 = np.zeros((nocc, nocc, nvir), dtype=ds_type) - vector1[p] = -1.0 - vector1 += np.einsum('ia,a->i', l1, t1[p, :]) - vector1 += 2 * np.einsum('ilcd,lcd->i', l2, t2[p, :, :, :]) - vector1 -= np.einsum('ilcd,ldc->i', l2, t2[p, :, :, :]) - - vector2[p, :, :] += -2. * l1 - vector2[:, p, :] += l1 - vector2 += 2 * np.einsum('c,ijcb->ijb', t1[p, :], l2) - vector2 -= np.einsum('c,jicb->ijb', t1[p, :], l2) - - return vector1, vector2 + result = np.zeros((nocc, nocc, nvir), dtype=ds_type) + result[p, :, :] += -2. * l1 + result[:, p, :] += l1 + result += 2 * np.einsum('c,ijcb->ijb', t1[p, :], l2) + result -= np.einsum('c,jicb->ijb', t1[p, :], l2) + return result else: - return -l1[:, p - nocc], -2 * l2[:, :, p - nocc, :] + l2[:, :, :, p - nocc] + return -2 * l2[:, :, p - nocc, :] + l2[:, :, :, p - nocc] def greens_e_vector_ip_rhf(cc, p): - return amplitudes_to_vector_ip(*greens_e_amplitudes_ip_rhf(cc.t1, cc.t2, cc.l1, cc.l2, p)) + return amplitudes_to_vector_ip( + greens_e_singles_ip_rhf(cc.t1, cc.t2, cc.l1, cc.l2, p), + greens_e_doubles_ip_rhf(cc.t1, cc.l1, cc.l2, p), + ) def greens_func_multiply(ham, vector, linear_part, **kwargs):