from mpi4py import MPI
import numpy as np
import scipy
from pyscf import lib, ao2mo
from pyscf.ao2mo import _ao2mo
from fcdmft.gw.mol.gw_exact import diagonalize_phrpa
from fcdmft.utils import cholesky
einsum = lib.einsum
rank = MPI.COMM_WORLD.Get_rank()
[docs]
def rpa_correlation_energy(mf, eri=None, Lpq=None, nw=40, exact=False):
"""RPA correlation energy
TODO: unrestricted.
Parameters
----------
mf : pyscf.scf.hf.RHF object
eri : double ndarray, optional
s1-symmetry two-electron integral in MO space, by default None
Lpq : double ndarray, optional
density fitting 3-center integral in MO basis, by default None
nw : int, optional
number of frequency point on imaginary axis, by default None
exact : bool, optional
use exact diagonalization, by default False
Returns
-------
ec : double
RPA correlation energy
"""
from fcdmft.rpa.mol.rpa import RPA, get_rpa_ecorr, _get_scaled_legendre_roots
assert eri is not None or Lpq is not None
def do_cd(full_eri):
try:
cd = cholesky.cholesky(full_eri, tau=1e-5, dimQ=250)
except:
if rank == 0:
print ('cderi tau=1e-5 fails, switch to tau=5e-5', flush=True)
cd = cholesky.cholesky(full_eri, tau=5e-5, dimQ=350)
cderi = cd.kernel()
return cderi
if Lpq is None:
nao = mf.mo_coeff.shape[0]
if len(eri.shape) != 4:
eri = ao2mo.restore(1, eri, nao)
cderi = do_cd(eri)
del eri
cderi = cderi.reshape(-1, nao, nao)
naux = cderi.shape[0]
if rank == 0:
print ('3-index ERI', cderi.shape)
Lpq_ao = cderi.copy()
del cderi
mo = np.asarray(mf.mo_coeff, order='F')
ijslice = (0, nao, 0, nao)
Lpq = _ao2mo.nr_e2(Lpq_ao, mo, ijslice, aosym='s1', mosym='s1', out=Lpq)
Lpq = Lpq.reshape(naux, nao, nao)
if exact is False:
rpa = RPA(mf)
freqs, wts = _get_scaled_legendre_roots(nw=nw)
if Lpq.shape[1] == mf.mo_coeff.shape[0]:
nocc = mf.mol.nelectron // 2
ec = get_rpa_ecorr(rpa=rpa, Lpq=Lpq[:,:nocc,nocc:], freqs=freqs, wts=wts)
else:
ec = get_rpa_ecorr(rpa=rpa, Lpq=Lpq, freqs=freqs, wts=wts)
else:
assert (Lpq.shape[1] == mf.mo_coeff.shape[0])
nocc = mf.mol.nelectron // 2
mo_energy = np.asarray(mf.mo_energy)
exci, _ = diagonalize_phrpa(nocc=nocc, mo_energy=mo_energy, Lpq=Lpq)
trace_A = 2.0 * einsum('Lia,Lia->', Lpq[:,:nocc,nocc:], Lpq[:,:nocc,nocc:])
trace_A += np.sum(np.asarray(mo_energy[None, nocc:] - mo_energy[:nocc, None]).reshape(-1))
sum_exci = np.sum(exci)
ec = 0.5 * (sum_exci - trace_A)
return ec
[docs]
def get_Lpq_emb_gamma(mf, lo_coeff, mf_cas):
"""Get Lpq 3-index integral for embedding problem from Gamma-point AO integral
Parameters
----------
mf : pyscf.scf.hf.RHF object
full mean-field object
lo_coeff : double ndarray
LO/EO coefficient, (nao_full, nao_cas)
mf_cas : pyscf.scf.hf.RHF object
embedding mean-field object
Returns
-------
Lpq : double ndarray
three-center density fitting matrix in mf_cas MO space, (naux_full, nmo_cas, nmo_cas)
"""
from fcdmft.rpa.pbc.rpa import RPA
myrpa = RPA(mf)
mo_coeff = np.dot(lo_coeff, mf_cas.mo_coeff)
Lpq = myrpa.ao2mo(mo_coeff=mo_coeff, fullmo=True)
return Lpq
[docs]
def diagonalize_pprpa_singlet(nocc, mo_energy, Lpq):
"""Diagonalize triplet ppRPA matrix.
Equation.21 in https://doi.org/10.1063/1.4828728.
Args:
nocc (int): number of occupied orbitals.
mo_energy (double array): orbital energy.
Lpq (double ndarray): three-center density fitting matrix in MO space.
Returns:
exci (double array): two-electron removal and addition energy.
vec (double ndarray): two-electron removal and addition eigenvector.
trace_A (double): trace of A matrix.
sum_exci (double): sum of two-electron addition energy.
"""
nmo = len(mo_energy)
nvir = nmo - nocc
mu = (mo_energy[nocc-1] + mo_energy[nocc]) * 0.5 # chemical potential
# get triangular part of full matrix, we need a>=b and i>=j
tri_row_o, tri_col_o = np.triu_indices(nocc)
tri_row_v, tri_col_v = np.triu_indices(nvir)
# A matrix: particle-particle block
A = einsum("Pac,Pbd->abcd", Lpq[:, nocc:, nocc:], Lpq[:, nocc:, nocc:])
A += einsum("Pad,Pbc->abcd", Lpq[:, nocc:, nocc:], Lpq[:, nocc:, nocc:])
A[np.diag_indices(nvir)] *= 1.0 / np.sqrt(2)
A = A.transpose(2, 3, 0, 1)
A[np.diag_indices(nvir)] *= 1.0 / np.sqrt(2)
A = A.transpose(2, 3, 0, 1)
A = A.reshape(nvir*nvir, nvir*nvir)
orb_sum = np.asarray(mo_energy[nocc:, None] + mo_energy[None, nocc:]).reshape(-1)
orb_sum -= 2.0 * mu
np.fill_diagonal(A, A.diagonal() + orb_sum)
A = A.reshape(nvir, nvir, nvir, nvir)
A = A[tri_row_v, tri_col_v, ...]
A = A[..., tri_row_v, tri_col_v]
# B matrix: particle-hole block
B = einsum("Pai,Pbj->abij", Lpq[:, nocc:, :nocc], Lpq[:, nocc:, :nocc])
B += einsum("Paj,Pbi->abij", Lpq[:, nocc:, :nocc], Lpq[:, nocc:, :nocc])
B[np.diag_indices(nvir)] *= 1.0 / np.sqrt(2)
B = B.transpose(2, 3, 0, 1)
B[np.diag_indices(nocc)] *= 1.0 / np.sqrt(2)
B = B.transpose(2, 3, 0, 1)
B = B[tri_row_v, tri_col_v, ...]
B = B[..., tri_row_o, tri_col_o]
# C matrix: hole-hole block
C = einsum("Pik,Pjl->ijkl", Lpq[:, :nocc, :nocc], Lpq[:, :nocc, :nocc])
C += einsum("Pil,Pjk->ijkl", Lpq[:, :nocc, :nocc], Lpq[:, :nocc, :nocc])
C[np.diag_indices(nocc)] *= 1.0 / np.sqrt(2)
C = C.transpose(2, 3, 0, 1)
C[np.diag_indices(nocc)] *= 1.0 / np.sqrt(2)
C = C.transpose(2, 3, 0, 1)
C = C.reshape(nocc*nocc, nocc*nocc)
orb_sum = np.asarray(mo_energy[:nocc, None] + mo_energy[None, :nocc]).reshape(-1)
orb_sum -= 2.0 * mu
np.fill_diagonal(C, C.diagonal() - orb_sum)
C = C.reshape(nocc, nocc, nocc, nocc)
C = C[tri_row_o, tri_col_o, ...]
C = C[..., tri_row_o, tri_col_o]
# combine A, B and C matrix as
# | A B |
# |B^T C|
M_upper = np.concatenate((A, B), axis=1)
M_lower = np.concatenate((B.T, C), axis=1)
M = np.concatenate((M_upper, M_lower), axis=0)
# M = MW, W is the metric matrix [[I,0],[0,-I]]
oo_dim = int((nocc + 1) * nocc / 2) # number of hole-hole pairs
vv_dim = int((nvir + 1) * nvir / 2) # number of particle-particle pairs
M[:][vv_dim:] *= -1.0
exci, vec = scipy.linalg.eigh(M)
trace_A = np.trace(A)
sum_exci = np.sum(exci[oo_dim:])
return exci, vec, trace_A, sum_exci
[docs]
def diagonalize_pprpa_triplet(nocc, mo_energy, Lpq):
"""Diagonalize triplet ppRPA matrix.
Equation.19 in https://doi.org/10.1063/1.4828728.
Args:
nocc (int): number of occupied orbitals.
mo_energy (double array): orbital energy.
Lpq (double ndarray): three-center density fitting matrix in MO space.
Returns:
exci (double array): two-electron removal and addition energy.
vec (double ndarray): two-electron removal and addition eigenvector.
trace_A (double): trace of A matrix.
sum_exci (double): sum of two-electron addition energy.
"""
nmo = len(mo_energy)
nvir = nmo - nocc
mu = (mo_energy[nocc-1] + mo_energy[nocc]) * 0.5 # chemical potential
# get triangular part (without diaognal) of full matrix, we need a>b and i>j
tri_row_o, tri_col_o = np.triu_indices(nocc, 1)
tri_row_v, tri_col_v = np.triu_indices(nvir, 1)
# A matrix: particle-particle block
A = einsum("Pac,Pbd->abcd", Lpq[:, nocc:, nocc:], Lpq[:, nocc:, nocc:])
A -= einsum("Pad,Pbc->abcd", Lpq[:, nocc:, nocc:], Lpq[:, nocc:, nocc:])
A = A.reshape(nvir*nvir, nvir*nvir)
orb_sum = np.asarray(mo_energy[nocc:, None] + mo_energy[None, nocc:]).reshape(-1)
orb_sum -= 2.0 * mu
np.fill_diagonal(A, A.diagonal() + orb_sum)
A = A.reshape(nvir, nvir, nvir, nvir)
A = A[tri_row_v, tri_col_v, ...]
A = A[..., tri_row_v, tri_col_v]
# B matrix: particle-hole block
B = einsum("Pai,Pbj->abij", Lpq[:, nocc:, :nocc], Lpq[:, nocc:, :nocc])
B -= einsum("Paj,Pbi->abij", Lpq[:, nocc:, :nocc], Lpq[:, nocc:, :nocc])
B = B[tri_row_v, tri_col_v, ...]
B = B[..., tri_row_o, tri_col_o]
# C matrix: hole-hole block
C = einsum("Pik,Pjl->ijkl", Lpq[:, :nocc, :nocc], Lpq[:, :nocc, :nocc])
C -= einsum("Pil,Pjk->ijkl", Lpq[:, :nocc, :nocc], Lpq[:, :nocc, :nocc])
C = C.reshape(nocc*nocc, nocc*nocc)
orb_sum = np.asarray(mo_energy[:nocc, None] + mo_energy[None, :nocc]).reshape(-1)
orb_sum -= 2.0 * mu
np.fill_diagonal(C, C.diagonal() - orb_sum)
C = C.reshape(nocc, nocc, nocc, nocc)
C = C[tri_row_o, tri_col_o, ...]
C = C[..., tri_row_o, tri_col_o]
# combine A, B and C matrix as
# | A B |
# |B^T C|
M_upper = np.concatenate((A, B), axis=1)
M_lower = np.concatenate((B.T, C), axis=1)
M = np.concatenate((M_upper, M_lower), axis=0)
# M = MW, W is the metric matrix [[I,0],[0,-I]]
oo_dim = int((nocc - 1) * nocc / 2) # number of hole-hole pairs
vv_dim = int((nvir - 1) * nvir / 2) # number of particle-particle pairs
M[:][vv_dim:] *= -1.0
exci, vec = scipy.linalg.eigh(M)
trace_A = np.trace(A)
sum_exci = np.sum(exci[oo_dim:])
return exci, vec, trace_A, sum_exci
[docs]
def pprpa_correlation_energy(mf, eri=None, Lpq=None):
"""Get ppRPA correlation energy.
Equation.4 in https://journals.aps.org/pra/abstract/10.1103/PhysRevA.88.030501.
Equation.24 in https://doi.org/10.1063/1.4828728.
TODO: unrestricted.
Args:
mf (pyscf.scf.hf.RHF): mean-field object.
eri (double ndarray): s1-symmetry two-electron integral in MO space.
Lpq (double ndarray): three-center density fitting matrix in MO space.
Returns:
ec (double): correlation energy.
"""
assert eri is not None or Lpq is not None
def do_cd(full_eri):
try:
cd = cholesky.cholesky(full_eri, tau=1e-5, dimQ=250)
except:
if rank == 0:
print ('cderi tau=1e-5 fails, switch to tau=5e-5', flush=True)
cd = cholesky.cholesky(full_eri, tau=5e-5, dimQ=350)
cderi = cd.kernel()
return cderi
if Lpq is None:
nao = mf.mo_coeff.shape[0]
if len(eri.shape) != 4:
eri = ao2mo.restore(1, eri, nao)
cderi = do_cd(eri)
del eri
cderi = cderi.reshape(-1, nao, nao)
naux = cderi.shape[0]
if rank == 0:
print ('3-index ERI', cderi.shape)
Lpq_ao = cderi.copy()
del cderi
mo = np.asarray(mf.mo_coeff, order='F')
ijslice = (0, nao, 0, nao)
Lpq = _ao2mo.nr_e2(Lpq_ao, mo, ijslice, aosym='s1', mosym='s1', out=Lpq)
Lpq = Lpq.reshape(naux, nao, nao)
assert (Lpq.shape[1] == mf.mo_coeff.shape[0])
nocc = mf.mol.nelectron // 2
mo_energy = np.asarray(mf.mo_energy)
trace_A_s, sum_exci_s = diagonalize_pprpa_singlet(nocc=nocc, mo_energy=mo_energy, Lpq=Lpq)[2:4]
trace_A_t, sum_exci_t = diagonalize_pprpa_triplet(nocc=nocc, mo_energy=mo_energy, Lpq=Lpq)[2:4]
trace_A = trace_A_s + trace_A_t * 3.0
sum_exci = sum_exci_s + sum_exci_t * 3.0
ec = sum_exci - trace_A
return ec