import time, sys, ctypes
import numpy as np
import numpy, math, scipy
from fcdmft.solver import gmres
from pyscf.lib import logger
from pyscf import lib, ao2mo, fci
from pyscf.fci.direct_spin1 import _unpack
libfci = lib.load_library('libfci')
[docs]
def greens_func_multiply_ip(ham, vector, linear_part, **kwargs):
return np.array(ham(vector.real, **kwargs) + 1j * ham(vector.imag, **kwargs) + linear_part * vector)
[docs]
def greens_func_multiply_ea(ham, vector, linear_part, **kwargs):
return np.array(-ham(vector.real, **kwargs) - 1j * ham(vector.imag, **kwargs) + linear_part * vector)
[docs]
def greens_func_multiply_ip_mor(ham, vector, linear_part):
return np.array(np.dot(ham, vector) + linear_part * vector)
[docs]
def greens_func_multiply_ea_mor(ham, vector, linear_part):
return np.array(-np.dot(ham, vector) + linear_part * vector)
[docs]
def contract_2e(fcivec, eri, norb, nelec, link_index=None):
fcivec = numpy.asarray(fcivec, order='C')
eri = ao2mo.restore(4, eri, norb)
link_indexa, link_indexb = _unpack(norb, nelec, link_index)
na, nlinka = link_indexa.shape[:2]
nb, nlinkb = link_indexb.shape[:2]
assert(fcivec.size == na*nb)
ci1 = numpy.empty_like(fcivec)
libfci.FCIcontract_2e_spin1(eri.ctypes.data_as(ctypes.c_void_p),
fcivec.ctypes.data_as(ctypes.c_void_p),
ci1.ctypes.data_as(ctypes.c_void_p),
ctypes.c_int(norb),
ctypes.c_int(na), ctypes.c_int(nb),
ctypes.c_int(nlinka), ctypes.c_int(nlinkb),
link_indexa.ctypes.data_as(ctypes.c_void_p),
link_indexb.ctypes.data_as(ctypes.c_void_p))
return ci1
[docs]
class FCIGF(object):
def __init__(self, myfci, mf, tol=1e-4, verbose=None):
self._fci = myfci
self._scf = mf
self.tol = tol
self.max_memory = mf.max_memory
if verbose:
self.verbose = verbose
else:
self.verbose = self._fci.verbose
self.stdout = sys.stdout
[docs]
def ipfci_mo(self, ps, qs, omega_list, broadening, MOR=False, omega_mor=None):
'''
Compute FCI IP Green's function in MO basis
Note: this is the fast routine without constructing the full FCI Hamiltonian
'''
mf = self._scf
myfci = self._fci
norb = len(mf.mo_energy)
nelec = mf.mol.nelec
neleca, nelecb = nelec
# 1e and 2e integrals in MO basis
h1e = mf.mo_coeff.T.dot(mf.get_hcore()).dot(mf.mo_coeff)
try:
eri = ao2mo.kernel(mf.mol, mf.mo_coeff)
except:
eri = ao2mo.kernel(mf._eri, mf.mo_coeff)
e_fci = myfci.energy(h1e, eri, myfci.ci, norb, nelec)
# make N-1 system h2e for calculating FCI H \dot c
h2e = myfci.absorb_h1e(h1e, eri, norb, (neleca-1,nelecb), 0.5)
# N-1 diagonal Hamiltonian for preconditioner
diag = self._fci.make_hdiag(h1e, eri, norb, (neleca-1,nelecb))
del eri
e_vector = list()
for q in qs:
e_vector.append(fci.addons.des_a(myfci.ci, norb, (neleca, nelecb), q).reshape(-1))
gfvals = np.zeros((len(ps), len(qs), len(omega_list)), dtype=np.complex128)
if MOR:
omega_comp = omega_mor
else:
omega_comp = omega_list
for ip, p in enumerate(ps):
b_vector = fci.addons.des_a(myfci.ci, norb, (neleca, nelecb), p).reshape(-1)
if MOR:
X_vec = np.zeros((len(b_vector), len(omega_mor)), dtype=np.complex128)
for iomega, curr_omega in enumerate(omega_comp):
def matr_multiply(vector, args=None):
return greens_func_multiply_ip(contract_2e, vector, curr_omega - e_fci - 1j * broadening,
eri=h2e, norb=norb, nelec=(neleca-1, nelecb))
diag_w = diag + curr_omega - e_fci - 1j*broadening
x0 = b_vector / diag_w
solver = gmres.GMRES(matr_multiply, b_vector, x0, diag_w, tol=self.tol)
cput1 = (time.process_time(), time.perf_counter())
sol = solver.solve().reshape(-1)
if MOR:
# Construct MOR subspace vectors
X_vec[:,iomega] = sol
cput1 = logger.timer(self, 'IPGF GMRES orbital p = %d/%d, freq w = %d/%d (%d iterations)'%(
ip+1,len(ps),iomega+1,len(omega_comp),solver.niter), *cput1)
x0 = sol
if not MOR:
for iq, q in enumerate(qs):
gfvals[ip,iq,iomega] = np.dot(e_vector[iq], sol)
if MOR:
# QR decomposition to orthorgonalize X_vec
S_vec, R_vec = scipy.linalg.qr(X_vec, mode='economic')
# Construct reduced order model
n_mor = len(omega_mor)
HS = np.zeros_like(S_vec)
for i in range(n_mor):
HS[:,i] = greens_func_multiply_ip(contract_2e, S_vec[:,i], 0.,
eri=h2e, norb=norb, nelec=(neleca-1, nelecb))
# Reduced Hamiltonian, b_vector, diag
H_mor = np.dot(S_vec.T.conj(), HS)
b_vector_mor = np.dot(S_vec.T.conj(), b_vector)
diag_mor = H_mor.diagonal()
for iomega in range(len(omega_list)):
curr_omega = omega_list[iomega]
def matr_multiply_mor(vector, args=None):
return greens_func_multiply_ip_mor(H_mor, vector, curr_omega - e_fci - 1j * broadening)
diag_w = diag_mor + curr_omega - e_fci - 1j*broadening
x0 = b_vector_mor / diag_w
solver = gmres.GMRES(matr_multiply_mor, b_vector_mor, x0, diag_w, tol=self.tol*0.01)
sol = solver.solve().reshape(-1)
x0 = sol
for iq, q in enumerate(qs):
e_vector_mor = np.dot(e_vector[iq], S_vec)
gfvals[ip,iq,iomega] = np.dot(e_vector_mor, sol)
return gfvals
[docs]
def eafci_mo(self, ps, qs, omega_list, broadening, MOR=False, omega_mor=None):
'''
Compute FCI EA Green's function in MO basis
Note: this is the fast routine without constructing the full FCI Hamiltonian
'''
mf = self._scf
myfci = self._fci
norb = len(mf.mo_energy)
nelec = mf.mol.nelec
neleca, nelecb = nelec
# 1e and 2e integrals in MO basis
h1e = mf.mo_coeff.T.dot(mf.get_hcore()).dot(mf.mo_coeff)
try:
eri = ao2mo.kernel(mf.mol, mf.mo_coeff)
except:
eri = ao2mo.kernel(mf._eri, mf.mo_coeff)
e_fci = myfci.energy(h1e, eri, myfci.ci, norb, nelec)
# make N+1 system h2e for calculating FCI H \dot c
h2e = myfci.absorb_h1e(h1e, eri, norb, (neleca+1,nelecb), 0.5)
# N+1 diagonal Hamiltonian for preconditioner
diag = self._fci.make_hdiag(h1e, eri, norb, (neleca+1,nelecb))
del eri
e_vector = list()
for q in qs:
e_vector.append(fci.addons.cre_a(myfci.ci, norb, (neleca, nelecb), q).reshape(-1))
gfvals = np.zeros((len(ps), len(qs), len(omega_list)), dtype=np.complex128)
if MOR:
omega_comp = omega_mor
else:
omega_comp = omega_list
for ip, p in enumerate(ps):
b_vector = fci.addons.cre_a(myfci.ci, norb, (neleca, nelecb), p).reshape(-1)
if MOR:
X_vec = np.zeros((len(b_vector), len(omega_mor)), dtype=np.complex128)
for iomega, curr_omega in enumerate(omega_comp):
def matr_multiply(vector, args=None):
return greens_func_multiply_ea(contract_2e, vector, curr_omega + e_fci + 1j * broadening,
eri=h2e, norb=norb, nelec=(neleca+1, nelecb))
diag_w = -diag + curr_omega + e_fci + 1j*broadening
x0 = b_vector / diag_w
solver = gmres.GMRES(matr_multiply, b_vector, x0, diag_w, tol=self.tol)
cput1 = (time.process_time(), time.perf_counter())
sol = solver.solve().reshape(-1)
if MOR:
# Construct MOR subspace vectors
X_vec[:,iomega] = sol
cput1 = logger.timer(self, 'EAGF GMRES orbital p = %d/%d, freq w = %d/%d (%d iterations)'%(
ip+1,len(ps),iomega+1,len(omega_comp),solver.niter), *cput1)
x0 = sol
if not MOR:
for iq, q in enumerate(qs):
gfvals[ip,iq,iomega] = np.dot(e_vector[iq], sol)
if MOR:
# QR decomposition to orthorgonalize X_vec
S_vec, R_vec = scipy.linalg.qr(X_vec, mode='economic')
# Construct reduced order model
n_mor = len(omega_mor)
HS = np.zeros_like(S_vec)
for i in range(n_mor):
HS[:,i] = -greens_func_multiply_ea(contract_2e, S_vec[:,i], 0.,
eri=h2e, norb=norb, nelec=(neleca+1, nelecb))
# Reduced Hamiltonian, b_vector, diag
H_mor = np.dot(S_vec.T.conj(), HS)
b_vector_mor = np.dot(S_vec.T.conj(), b_vector)
diag_mor = H_mor.diagonal()
for iomega in range(len(omega_list)):
curr_omega = omega_list[iomega]
def matr_multiply_mor(vector, args=None):
return greens_func_multiply_ea_mor(H_mor, vector, curr_omega + e_fci + 1j * broadening)
diag_w = -diag_mor + curr_omega + e_fci + 1j*broadening
x0 = b_vector_mor / diag_w
solver = gmres.GMRES(matr_multiply_mor, b_vector_mor, x0, diag_w, tol=self.tol*0.01)
sol = solver.solve().reshape(-1)
x0 = sol
for iq, q in enumerate(qs):
e_vector_mor = np.dot(e_vector[iq], S_vec)
gfvals[ip,iq,iomega] = np.dot(e_vector_mor, sol)
return gfvals
[docs]
def ipfci_mo_slow(self, ps, qs, omega_list, broadening, np_max=1e5):
'''
Slow version of computing FCI IP Green's function in MO basis
Warning: this routine constructs FCI Hamiltonian matrix and only works for small active space
'''
mf = self._scf
myfci = self._fci
norb = len(mf.mo_energy)
nelec = mf.mol.nelec
neleca, nelecb = nelec
# 1e and 2e integrals in MO basis
h1e = mf.mo_coeff.T.dot(mf.get_hcore()).dot(mf.mo_coeff)
try:
eri = ao2mo.kernel(mf.mol, mf.mo_coeff)
except:
eri = ao2mo.kernel(mf._eri, mf.mo_coeff)
# FCI ground-state electronic energy
ene_fci = myfci.energy(h1e, eri, myfci.ci, norb, nelec)
vec_fci = myfci.ci.reshape(-1)
size = vec_fci.size
nelec_ip = (nelec[0] - 1, nelec[1])
nelec_ea = (nelec[0] + 1, nelec[1])
omegas = numpy.asarray(omega_list)
nomega = len(omegas)
ps = numpy.asarray(ps)
qs = numpy.asarray(qs)
np = len(ps)
nq = len(qs)
# Set up the IP manybody Hamiltonian
link_index_ip = _unpack(norb, nelec_ip, None, spin=None)
na_ip = link_index_ip[0].shape[0]
nb_ip = link_index_ip[1].shape[0]
size_ip = na_ip * nb_ip
hdiag_ip = myfci.make_hdiag(h1e, eri, norb, nelec_ip)
assert hdiag_ip.shape == (size_ip, )
# Build the RHS and LHS of the response equation
bps = numpy.asarray([fci.addons.des_a(vec_fci, norb, nelec, p).reshape(-1) for p in ps]).reshape(np, size_ip)
eqs = numpy.asarray([fci.addons.des_a(vec_fci, norb, nelec, q).reshape(-1) for q in qs]).reshape(nq, size_ip)
# h2e_ip = myfci.absorb_h1e(h1e, eri, norb, nelec_ip, fac=0.5)
if size_ip * size_ip * 8 / 1024**3 > self.max_memory:
raise ValueError("Not enough memory for FCI IP Hamiltonian.")
h_ip = fci.direct_spin1.pspace(h1e, eri, norb, nelec_ip, hdiag=hdiag_ip, np=size_ip)[1]
assert h_ip.shape == (size_ip, size_ip)
def gen_gfn(omega):
omega_e0_eta_ip = omega - ene_fci - 1j * eta
h_ip_omega = h_ip + omega_e0_eta_ip * numpy.eye(size_ip)
# h_ip_omega_inv = numpy.linalg.inv(h_ip_omega)
# xps = numpy.dot(h_ip_omega_inv, bps.T)
xps = numpy.linalg.solve(h_ip_omega, bps.T).T
return numpy.dot(xps, eqs.T)
cput0 = (lib.logger.process_clock(), lib.logger.perf_counter())
gfns_ip = numpy.asarray([gen_gfn(omega) for omega in omegas]).reshape(nomega, np, nq)
lib.logger.timer(self, 'Solving GF-IP', *cput0)
return gfns_ip.transpose(1, 2, 0)
[docs]
def eafci_mo_slow(self, ps, qs, omega_list, broadening, np_max=1e5):
'''
Slow version of computing FCI EA Green's function in MO basis
Warning: this routine constructs FCI Hamiltonian matrix and only works for small active space
'''
mf = self._scf
myfci = self._fci
norb = len(mf.mo_energy)
nelec = mf.mol.nelec
neleca, nelecb = nelec
# 1e and 2e integrals in MO basis
h1e = mf.mo_coeff.T.dot(mf.get_hcore()).dot(mf.mo_coeff)
try:
eri = ao2mo.kernel(mf.mol, mf.mo_coeff)
except:
eri = ao2mo.kernel(mf._eri, mf.mo_coeff)
# FCI ground-state electronic energy
ene_fci = myfci.energy(h1e, eri, myfci.ci, norb, nelec)
vec_fci = myfci.ci.reshape(-1)
size = vec_fci.size
nelec_ip = (nelec[0] - 1, nelec[1])
nelec_ea = (nelec[0] + 1, nelec[1])
omegas = numpy.asarray(omega_list)
nomega = len(omegas)
ps = numpy.asarray(ps)
qs = numpy.asarray(qs)
np = len(ps)
nq = len(qs)
# Set up the EA manybody Hamiltonian
link_index_ea = _unpack(norb, nelec_ea, None, spin=None)
na_ea = link_index_ea[0].shape[0]
nb_ea = link_index_ea[1].shape[0]
size_ea = na_ea * nb_ea
hdiag_ea = myfci.make_hdiag(h1e, eri, norb, nelec_ea)
assert hdiag_ea.shape == (size_ea, )
bps = numpy.asarray([fci.addons.cre_a(vec_fci, norb, nelec, p).reshape(-1) for p in ps]).reshape(np, size_ea)
eqs = numpy.asarray([fci.addons.cre_a(vec_fci, norb, nelec, q).reshape(-1) for q in qs]).reshape(nq, size_ea)
# h2e_ea = myfci.absorb_h1e(h1e, eri, norb, nelec_ea, fac=0.5)
if size_ea * size_ea * 8 / 1024**3 > self.max_memory:
raise ValueError("Not enough memory for FCI ea Hamiltonian.")
h_ea = fci.direct_spin1.pspace(h1e, eri, norb, nelec_ea, hdiag=hdiag_ea, np=size_ea)[1]
assert h_ea.shape == (size_ea, size_ea)
def gen_gfn(omega):
omega_e0_eta_ea = omega + ene_fci + 1j * eta
h_ea_omega = - h_ea + omega_e0_eta_ea * numpy.eye(size_ea)
# h_ea_omega_inv = numpy.linalg.inv(h_ea_omega)
# xps = numpy.dot(h_ea_omega_inv, bps.T)
xps = numpy.linalg.solve(h_ea_omega, bps.T).T
return numpy.dot(xps, eqs.T)
gfns_ea = numpy.asarray([gen_gfn(omega) for omega in omegas]).reshape(nomega, np, nq)
return gfns_ea.transpose(1, 2, 0)
if __name__ == '__main__':
from pyscf import gto, scf
mol = gto.M(
atom = 'H 0 0 0; Li 0 0 1.1',
basis = 'sto-3g',
verbose = 5,
)
mf = scf.RHF(mol)
mf.kernel()
myfci = fci.FCI(mf)
e, c = myfci.kernel()
myfcigf = FCIGF(myfci, mf, tol=1e-8)
omega_list = np.linspace(-0.5, 0.5, 201)
eta = 0.01
nmo = len(mf.mo_energy)
ps = range(nmo)
qs = range(nmo)
gf_fci_ip_1 = myfcigf.ipfci_mo_slow(ps, qs, omega_list, eta).conj()
gf_fci_ea_1 = myfcigf.eafci_mo_slow(ps, qs, omega_list, eta)
gf_fci_ip_2 = myfcigf.ipfci_mo(ps, qs, omega_list, eta).conj()
gf_fci_ea_2 = myfcigf.eafci_mo(ps, qs, omega_list, eta)
gf_fci = gf_fci_ip_2 + gf_fci_ea_2
assert np.linalg.norm(gf_fci_ip_1 - gf_fci_ip_2) < 1e-4
assert np.linalg.norm(gf_fci_ea_1 - gf_fci_ea_2) < 1e-4
print ('--- FCI GF ---')
for iomega, omega in enumerate(omega_list):
print("% 8.4f, %8.4f"%(omega, -np.trace(gf_fci[:,:,iomega].imag)/np.pi))
omega_mor_ip = np.linspace(-0.5, 0., 20)
omega_mor_ea = np.linspace(0., 0.5, 20)
gf_fci_ip = myfcigf.ipfci_mo(ps, qs, omega_list, eta, MOR=True, omega_mor=omega_mor_ip).conj()
gf_fci_ea = myfcigf.eafci_mo(ps, qs, omega_list, eta, MOR=True, omega_mor=omega_mor_ea)
gf_fci = gf_fci_ip + gf_fci_ea
print ('--- FCI MOR GF ---')
for iomega, omega in enumerate(omega_list):
print("% 8.4f, %8.4f"%(omega, -np.trace(gf_fci[:,:,iomega].imag)/np.pi))