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
import itertools
from mpi4py import MPI
import math
from functools import lru_cache
libfci = lib.load_library('libfci')
'''
MPI FCI Green's function
'''
rank = MPI.COMM_WORLD.Get_rank()
size = MPI.COMM_WORLD.Get_size()
comm = MPI.COMM_WORLD
[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
MPI parallelize over orbitals
'''
if MOR:
return self.ipfci_mo_mor(ps, qs, omega_list, broadening, omega_mor)
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))
omega_comp = omega_list
pws = list(itertools.product(ps, omega_comp))
global_p_omegas = list(zip(range(len(pws)), pws))
WSZ = len(global_p_omegas)
segsize = math.ceil(len(global_p_omegas) / size)
local_p_omegas = global_p_omegas[rank * segsize : min((rank + 1) * segsize, WSZ)]
gfvals = np.zeros((segsize, len(qs)), dtype=np.complex128)
for ipw, (gipw, (p, curr_omega)) in enumerate(local_p_omegas):
#p = ps[ip]
@lru_cache(maxsize=1)
def take_from_p(p):
return fci.addons.des_a(myfci.ci, norb, (neleca, nelecb), p).reshape(-1)
b_vector = take_from_p(p)
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)
cput1 = logger.timer(self, 'IPGF GMRES orbital ((p,w) = %d / %d) (%d iterations)'%(
gipw+1,WSZ,solver.niter), *cput1)
x0 = sol
for iq, q in enumerate(qs):
gfvals[ipw,iq] = np.dot(e_vector[iq], sol)
gfvals_all = np.empty((size * segsize, len(qs)), dtype=np.complex128)
comm.Barrier()
comm.Allgather([gfvals, MPI.DOUBLE_COMPLEX], [gfvals_all, MPI.DOUBLE_COMPLEX])
gfvals_all.resize((WSZ, len(qs)))
gfvals_all = np.reshape(gfvals_all, (len(ps), len(omega_list), len(qs)))
gfvals_all = np.transpose(gfvals_all, (0, 2, 1)).copy()
return gfvals_all
[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
'''
if MOR:
return self.eafci_mo_mor(ps, qs, omega_list, broadening, omega_mor)
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))
comm.Barrier()
omega_comp = omega_list
pws = list(itertools.product(ps, omega_comp))
global_p_omegas = list(zip(range(len(pws)), pws))
WSZ = len(global_p_omegas)
segsize = math.ceil(len(global_p_omegas) / size)
local_p_omegas = global_p_omegas[rank * segsize : min((rank + 1) * segsize, WSZ)]
gfvals = np.zeros((segsize, len(qs)), dtype=np.complex128)
for ipw, (gipw, (p, curr_omega)) in enumerate(local_p_omegas):
#p = ps[ip]
@lru_cache(maxsize=1)
def add_to_p(p):
return fci.addons.cre_a(myfci.ci, norb, (neleca, nelecb), p).reshape(-1)
b_vector = add_to_p(p)
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)
cput1 = logger.timer(self, 'EAGF GMRES orbital ((p,w) = %d / %d) (%d iterations)'%(
gipw+1,WSZ,solver.niter), *cput1)
x0 = sol
for iq, q in enumerate(qs):
gfvals[ipw,iq] = np.dot(e_vector[iq], sol)
gfvals_all = np.empty((size * segsize, len(qs)), dtype=np.complex128)
comm.Barrier()
comm.Allgather([gfvals, MPI.DOUBLE_COMPLEX], [gfvals_all, MPI.DOUBLE_COMPLEX])
gfvals_all.resize((WSZ, len(qs)))
gfvals_all = np.reshape(gfvals_all, (len(ps), len(omega_list), len(qs)))
gfvals_all = np.transpose(gfvals_all, (0, 2, 1)).copy()
return gfvals_all
[docs]
def ipfci_mo_mor(self, ps, qs, omega_list, broadening, omega_mor):
'''
Compute FCI IP Green's function in MO basis
Note: this is the fast routine without constructing the full FCI Hamiltonian
MPI parallelize over orbitals
'''
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 = [fci.addons.des_a(myfci.ci, norb, (neleca, nelecb), q).reshape(-1)
for q in qs]
comm.Barrier()
segsize = math.ceil(len(ps) / size)
local_ps = ps[rank * segsize : min((rank + 1) * segsize, len(ps))]
gfvals = np.zeros((segsize, len(qs), len(omega_list)), dtype=np.complex128)
for ip, p in enumerate(local_ps):
b_vector = fci.addons.des_a(myfci.ci, norb, (neleca, nelecb), p).reshape(-1)
X_vec = np.zeros((len(b_vector), len(omega_mor)), dtype=np.complex128)
for iomega, curr_omega in enumerate(omega_mor):
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)
# 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_mor),solver.niter), *cput1)
x0 = sol
# 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)
gfvals_all = np.empty((size * segsize, len(qs), len(omega_list)), dtype=np.complex128)
comm.Barrier()
comm.Allgather([gfvals, MPI.DOUBLE_COMPLEX], [gfvals_all, MPI.DOUBLE_COMPLEX])
gfvals_all.resize((len(ps), len(qs), len(omega_list)))
return gfvals_all
[docs]
def eafci_mo_mor(self, ps, qs, omega_list, broadening, omega_mor):
'''
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 = [fci.addons.cre_a(myfci.ci, norb, (neleca, nelecb), q).reshape(-1)
for q in qs]
comm.Barrier()
segsize = math.ceil(len(ps) / size)
local_ps = ps[rank * segsize : min((rank + 1) * segsize, len(ps))]
gfvals = np.zeros((segsize, len(qs), len(omega_list)), dtype=np.complex128)
for ip, p in enumerate(local_ps):
b_vector = fci.addons.cre_a(myfci.ci, norb, (neleca, nelecb), p).reshape(-1)
X_vec = np.zeros((len(b_vector), len(omega_mor)), dtype=np.complex128)
for iomega, curr_omega in enumerate(omega_mor):
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 True:
# 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_mor),solver.niter), *cput1)
x0 = sol
# 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)
gfvals_all = np.empty((size * segsize, len(qs), len(omega_list)), dtype=np.complex128)
comm.Barrier()
comm.Allgather([gfvals, MPI.DOUBLE_COMPLEX], [gfvals_all, MPI.DOUBLE_COMPLEX])
gfvals_all.resize((len(ps), len(qs), len(omega_list)))
return gfvals_all