#!/usr/bin/env python
# Copyright 2014-2021 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Author: Tianyu Zhu <zhutianyu1991@gmail.com>
#
"""
Periodic k-point spin-unrestricted random phase approximation
(direct RPA/dRPA in chemistry) with N^4 scaling
Method:
Main routines are based on GW-AC method described in:
T. Zhu and G.K.-L. Chan, J. Chem. Theory. Comput. 17, 727-741 (2021)
X. Ren et al., New J. Phys. 14, 053020 (2012)
"""
import time, os
import numpy as np
import scipy.linalg.blas as blas
from pyscf import lib
from pyscf.lib import logger, temporary_env
from pyscf.ao2mo import _ao2mo
from pyscf.ao2mo.incore import _conc_mos
from pyscf.pbc import df, dft, scf
from pyscf.pbc.cc.kccsd_uhf import get_nocc, get_nmo, get_frozen_mask
from fcdmft.gw.minimax.minimax_rpa import get_frequency_minimax_grid
from fcdmft.ac.grids import _get_scaled_legendre_roots
from fcdmft.gw.pbc.kugw_ac import get_rho_response, get_rho_response_metal, get_rho_response_head, \
get_rho_response_wing, get_qij
from fcdmft.rpa.pbc.krpa import KRPA, rho_accum_inner, get_rpa_ecorr_w
from fcdmft.rpa.mol.urpa import get_idx_metal
from fcdmft.utils.kpts import get_kconserv_ria_efficient
from mpi4py import MPI
rank = MPI.COMM_WORLD.Get_rank()
size = MPI.COMM_WORLD.Get_size()
comm = MPI.COMM_WORLD
[docs]
def kernel(rpa, mo_energy, mo_coeff, nw=None):
"""RPA correlation and total energy
Parameters
----------
rpa : KURPA
rpa object
mo_energy : double array
molecular orbital energies
mo_coeff : double ndarray
molecular orbital coefficients
nw : int, optional
number of frequency point on imaginary axis, by default None
Returns
-------
e_tot : float
RPA total energy
e_hf : float
HF energy (exact exchange for given mo_coeff)
e_corr : float
RPA correlation energy
"""
assert rpa.frozen == 0 or rpa.frozen is None
# Compute HF exchange energy (EXX)
dm = rpa._scf.make_rdm1()
uhf = scf.KUHF(rpa.mol, rpa.kpts, exxdiv=rpa._scf.exxdiv)
uhf.verbose = 0
if hasattr(rpa._scf, 'sigma'):
uhf = scf.addons.smearing_(uhf, sigma=rpa._scf.sigma, method=rpa._scf.smearing_method)
uhf.with_df = rpa._scf.with_df
with temporary_env(rpa.with_df, verbose=0), temporary_env(rpa.mol, verbose=0):
e_hf = uhf.energy_elec(dm)[0]
e_hf += rpa._scf.energy_nuc()
is_metal = hasattr(rpa._scf, 'sigma')
# Turn off FC for metals and outcore
if is_metal and rpa.fc:
if rank == 0:
logger.warn(rpa, 'FC not available for metals - setting rpa.fc to False')
rpa.fc = False
if rpa.fc and rpa.outcore:
if rank == 0:
logger.warn(rpa, 'FC not available for outcore - setting rpa.fc to False')
rpa.fc = False
# Grids for integration on imaginary axis
freqs, wts = rpa.get_grids(nw=nw, mo_energy=mo_energy)
# Compute RPA correlation energy
if rpa.outcore:
if is_metal:
e_corr = get_rpa_ecorr_outcore_metal(rpa, freqs, wts)
else:
e_corr = get_rpa_ecorr_outcore(rpa, freqs, wts)
else:
e_corr = get_rpa_ecorr(rpa, freqs, wts)
# Compute total energy
e_tot = e_hf + e_corr
if rank == 0 and rpa.verbose >= logger.DEBUG:
logger.debug(rpa, ' RPA total energy = %s', e_tot)
logger.debug(rpa, ' EXX energy = %s, RPA corr energy = %s', e_hf, e_corr)
return e_tot, e_hf, e_corr
[docs]
def get_rpa_ecorr(rpa, freqs, wts):
"""Compute RPA correlation energy.
Parameters
----------
rpa : KURPA
rpa object
freqs : double 1d array
frequency grid
wts : double 1d array
weight of grids
Returns
-------
e_corr : double
correlation energy
"""
mo_energy = np.array(rpa._scf.mo_energy)
mo_coeff = np.array(rpa._scf.mo_coeff)
nmoa, nmob = rpa.nmo
nkpts = rpa.nkpts
kpts = rpa.kpts
nw = len(freqs)
mydf = rpa.with_df
mo_occ = rpa.mo_occ
# possible kpts shift
kscaled = rpa.mol.get_scaled_kpts(kpts)
kscaled -= kscaled[0]
is_metal = hasattr(rpa._scf, 'sigma')
segsize = nkpts // size
if rank >= size - (nkpts - segsize * size):
start = rank * segsize + rank - (size - (nkpts - segsize * size))
stop = min(nkpts, start + segsize + 1)
else:
start = rank * segsize
stop = min(nkpts, start + segsize)
if rpa.fc:
qij_a, qij_b, q_abs, nq_pts = rpa.get_q_mesh(mo_energy, mo_coeff)
e_corr = 0j
# Precompute k-conservation table
# Given k-point indices (kL, i), kconserv_table[kshift,i] contains
# the index j that satisfies momentum conservation,
# (k(i) - k(j) - k(kL)) \dot a = 2n\pi
# i.e.
# - ki + kj + kL = G
kconserv_table = get_kconserv_ria_efficient(rpa.mol, kpts)
cderiarr = mydf.cderi_array()
for kL in range(start, stop):
# Lij: (2, ki, L, i, j) for looping every kL
# Lij = np.zeros((2,nkpts,naux,nmoa,nmoa),dtype=np.complex128)
Lij = []
# kidx: save kj that conserves with kL and ki (-ki+kj+kL=G)
# kidx_r: save ki that conserves with kL and kj (-ki+kj+kL=G)
kidx = np.zeros((nkpts), dtype=np.int64)
kidx_r = np.zeros((nkpts), dtype=np.int64)
for i, kpti in enumerate(kpts):
j = kconserv_table[kL, i]
kptj = kpts[j]
# Find (ki,kj) that satisfies momentum conservation with kL
kconserv = -kscaled[i] + kscaled[j] + kscaled[kL]
assert np.linalg.norm(np.round(kconserv) - kconserv) < 1e-12 # kidx[i] = j
kidx[i] = j
kidx_r[j] = i
logger.debug(rpa, 'Read Lpq (kL: %s / %s, ki: %s, kj: %s @ Rank %d)' % (kL + 1, nkpts, i, j, rank))
Lij_out_a = None
Lij_out_b = None
# Read (L|pq) and ao2mo transform to (L|ij)
Lpq = cderiarr.load(kpti, kptj)
if Lpq.shape[-1] == (nmoa * (nmoa + 1)) // 2:
Lpq = lib.unpack_tril(Lpq).reshape(-1, nmoa**2)
else:
Lpq = Lpq.reshape(-1, nmoa**2)
Lpq = Lpq.astype(np.complex128)
moija, ijslicea = _conc_mos(mo_coeff[0, i], mo_coeff[0, j])[2:]
moijb, ijsliceb = _conc_mos(mo_coeff[1, i], mo_coeff[1, j])[2:]
Lij_out_a = _ao2mo.r_e2(Lpq, moija, ijslicea, tao=[], ao_loc=None, out=Lij_out_a)
Lij_out_b = _ao2mo.r_e2(Lpq, moijb, ijsliceb, tao=[], ao_loc=None, out=Lij_out_b)
Lij.append(np.asarray((Lij_out_a.reshape(-1, nmoa, nmoa), Lij_out_b.reshape(-1, nmob, nmob))))
Lij = np.asarray(Lij)
naux = Lij.shape[2]
if is_metal is False:
Lia = [
np.ascontiguousarray(Lij[:, 0, :, : rpa.nocc[0], rpa.nocc[0] :]),
np.ascontiguousarray(Lij[:, 1, :, : rpa.nocc[1], rpa.nocc[1] :]),
]
for w in range(nw):
# body polarizability
if is_metal:
Pi = get_rho_response_metal(freqs[w], mo_energy, mo_occ, Lij, kidx)
else:
Pi = get_rho_response(freqs[w], rpa.nocc, mo_energy, Lia, kidx)
if kL == 0 and rpa.fc:
for iq in range(nq_pts):
# head Pi_00
Pi_00 = get_rho_response_head(freqs[w], mo_energy, (qij_a[iq], qij_b[iq]))
Pi_00 = 4.0 * np.pi / np.linalg.norm(q_abs[iq]) ** 2 * Pi_00
# wings Pi_P0
Pi_P0 = get_rho_response_wing(freqs[w], mo_energy, Lia, (qij_a[iq], qij_b[iq]))
Pi_P0 = np.sqrt(4.0 * np.pi) / np.linalg.norm(q_abs[iq]) * Pi_P0
# assemble Pi
Pi_fc = np.zeros((naux + 1, naux + 1), dtype=Pi.dtype)
Pi_fc[0, 0] = Pi_00
Pi_fc[0, 1:] = Pi_P0.conj()
Pi_fc[1:, 0] = Pi_P0
Pi_fc[1:, 1:] = Pi
# First, compute ec_w = Tr(Pi) + |log(det(I-Pi))|
ec_w = np.trace(Pi_fc)
# The following two lines are equivalent to
# Pi = np.eye(naux) - Pi
blas.zdscal(-1.0, Pi_fc.ravel(), overwrite_x=1)
np.fill_diagonal(Pi_fc, np.diagonal(Pi_fc) + 1.0)
ec_w += np.linalg.slogdet((Pi_fc))[1]
e_corr += 1.0 / (2.0 * np.pi) * 1.0 / nkpts * 1.0 / nq_pts * ec_w * wts[w]
else:
# First, compute ec_w = Tr(Pi) + |log(det(I-Pi))|
ec_w = np.trace(Pi)
# The following two lines are equivalent to
# Pi = np.eye(naux) - Pi
blas.zdscal(-1.0, Pi.ravel(), overwrite_x=1)
np.fill_diagonal(Pi, np.diagonal(Pi) + 1.0)
ec_w += np.linalg.slogdet((Pi))[1]
e_corr += 1.0 / (2.0 * np.pi) * 1.0 / nkpts * ec_w * wts[w]
comm.Barrier()
ecorr_gather = comm.gather(e_corr)
if rank == 0:
e_corr = np.sum(ecorr_gather)
comm.Barrier()
e_corr = comm.bcast(e_corr, root=0)
return e_corr.real
[docs]
def get_rpa_ecorr_outcore(rpa, freqs, wts):
"""Low-memory routine to compute RPA correlation energy.
Parameters
----------
rpa : KURPA
rpa object
freqs : double 1d array
frequency grid
wts : double 1d array
weight of grids
Returns
-------
e_corr : double
correlation energy
"""
mo_energy = np.array(rpa._scf.mo_energy)
mo_coeff = np.array(rpa._scf.mo_coeff)
nmoa = rpa.nmo[0]
nkpts = rpa.nkpts
kpts = rpa.kpts
nw = len(freqs)
mydf = rpa.with_df
# possible kpts shift
kscaled = rpa.mol.get_scaled_kpts(kpts)
kscaled -= kscaled[0]
segsize = nkpts // size
if rank >= size - (nkpts - segsize * size):
start = rank * segsize + rank - (size - (nkpts - segsize * size))
stop = min(nkpts, start + segsize + 1)
else:
start = rank * segsize
stop = min(nkpts, start + segsize)
if rpa.fc:
qij_a, qij_b, q_abs, nq_pts = rpa.get_q_mesh(mo_energy, mo_coeff)
e_corr = 0j
# Precompute k-conservation table
# Given k-point indices (kL, i), kconserv_table[kshift,i] contains
# the index j that satisfies momentum conservation,
# (k(i) - k(j) - k(kL)) \dot a = 2n\pi
# i.e.
# - ki + kj + kL = G
kconserv_table = get_kconserv_ria_efficient(rpa.mol, kpts)
cderiarr = mydf.cderi_array()
for kL in range(start, stop):
Pi = None
kidx = np.zeros((nkpts), dtype=np.int64)
kidx_r = np.zeros((nkpts), dtype=np.int64)
for s in range(2):
nseg = rpa.nocc[s] // rpa.segsize + 1
for iseg in range(nseg):
orb_start = iseg * rpa.segsize
orb_end = min((iseg + 1) * rpa.segsize, rpa.nocc[s])
if orb_end == orb_start:
continue
norb_this_iter = orb_end - orb_start
for i, kpti in enumerate(kpts):
j = kconserv_table[kL, i]
kptj = kpts[j]
# Find (ki,kj) that satisfies momentum conservation with kL
kconserv = -kscaled[i] + kscaled[j] + kscaled[kL]
assert np.linalg.norm(np.round(kconserv) - kconserv) < 1e-12 # kidx[i] = j
kidx[i] = j
kidx_r[j] = i
logger.debug(rpa, 'Read Lpq (kL: %s / %s, ki: %s, kj: %s @ Rank %d)' % (kL + 1, nkpts, i, j, rank))
# Read (L|pq) and ao2mo transform to (L|ij)
Lpq = cderiarr.load(kpti, kptj)
if Lpq.shape[-1] == (nmoa * (nmoa + 1)) // 2:
Lpq = lib.unpack_tril(Lpq).reshape(-1, nmoa**2)
else:
Lpq = Lpq.reshape(-1, nmoa**2)
Lpq = Lpq.astype(np.complex128)
naux = Lpq.shape[0]
moij, ijslice = _conc_mos(mo_coeff[s, i], mo_coeff[s, j])[2:]
ijslice = (orb_start, orb_end, rpa.nmo[s] + rpa.nocc[s], 2 * rpa.nmo[s])
Lij_slice = _ao2mo.r_e2(Lpq, moij, ijslice, tao=[], ao_loc=None)
Lij_slice = Lij_slice.reshape(naux, norb_this_iter, rpa.nmo[s] - rpa.nocc[s])
if Pi is None:
Pi = np.zeros((nw, naux, naux), dtype=np.complex128)
eia = mo_energy[s, i][orb_start:orb_end, None] - mo_energy[s, j][None, rpa.nocc[s] :]
for w in range(nw):
rho_accum_inner(Pi[w], eia, freqs[w], Lij_slice, alpha=2.0 / nkpts)
for w in range(nw):
e_corr += get_rpa_ecorr_w(Pi[w], wts[w])
e_corr = e_corr.real
e_corr = comm.allreduce(e_corr, op=MPI.SUM)
e_corr *= 1.0 / (2.0 * np.pi) / nkpts
return e_corr
[docs]
class KURPA(KRPA):
[docs]
def dump_flags(self, verbose=None):
log = logger.Logger(self.stdout, self.verbose)
log.info('')
log.info('******** %s ********', self.__class__)
log.info('method = %s', self.__class__.__name__)
nocca, noccb = self.nocc
nmoa, nmob = self.nmo
nvira = nmoa - nocca
nvirb = nmob - noccb
nkpts = self.nkpts
log.info(
'RPA (nocca, noccb) = (%d, %d), (nvira, nvirb) = (%d, %d), nkpts = %d', nocca, noccb, nvira, nvirb, nkpts
)
if self.frozen is not None and self.frozen > 0:
log.info('frozen orbitals %s', str(self.frozen))
log.info('grid type = %s', self.grids_alg)
log.info('outcore mode = %s', self.outcore)
if self.outcore is True:
log.info('outcore segment size = %d', self.segsize)
log.info('RPA finite size corrections = %s', self.fc)
log.info('')
return
@property
def nocc(self):
mo_occ = self._scf.mo_occ
return (int(np.sum(mo_occ[0][0])), int(np.sum(mo_occ[1][0])))
@nocc.setter
def nocc(self, n):
self._nocc = n
@property
def nmo(self):
return (len(self._scf.mo_energy[0][0]), len(self._scf.mo_energy[1][0]))
@nmo.setter
def nmo(self, n):
self._nmo = n
get_nocc = get_nocc
get_nmo = get_nmo
get_frozen_mask = get_frozen_mask
[docs]
def kernel(self, mo_energy=None, mo_coeff=None, nw=None):
"""RPA correlation and total energy
Calculated total energy, HF energy and RPA correlation energy
are stored in self.e_tot, self.e_hf, self.e_corr
Parameters
----------
mo_energy : double array
molecular orbital energies
mo_coeff : double ndarray
molecular orbital coefficients
nw : int, optional
number of frequency point on imaginary axis, by default None
Returns
-------
e_tot : float
RPA total energy
e_hf : float
HF energy (exact exchange for given mo_coeff)
e_corr : float
RPA correlation energy
"""
if mo_coeff is None:
mo_coeff = np.array(self._scf.mo_coeff)
if mo_energy is None:
mo_energy = np.array(self._scf.mo_energy)
nmoa = self.nmo[0]
naux = self.with_df.get_naoaux()
nkpts = self.nkpts
mem_incore = (3 * nkpts * nmoa**2 * naux) * 16 / 1e6
mem_now = lib.current_memory()[0]
if mem_incore + mem_now > 0.99 * self.max_memory:
if rank == 0:
logger.warn(self, 'Memory may not be enough!')
raise NotImplementedError
cput0 = (time.process_time(), time.perf_counter())
if rank == 0:
self.dump_flags()
self.e_tot, self.e_hf, self.e_corr = kernel(self, mo_energy, mo_coeff, nw=nw)
if rank == 0:
logger.timer(self, 'RPA', *cput0)
return self.e_tot, self.e_hf, self.e_corr
[docs]
def get_grids(self, alg=None, nw=None, mo_energy=None):
"""Generate grids for integration.
Parameters
----------
alg : str, optional
algorithm for generating grids, by default None
nw : int, optional
number of grids, by default None
mo_energy : double 3d array, optional
orbital energy, used for minimax grids, by default None
Returns
-------
freqs : double 1d array
frequency grid
wts : double 1d array
weight of grids
"""
if alg is None:
alg = self.grids_alg
if mo_energy is None:
mo_energy = np.array(self._scf.mo_energy)
if alg == 'legendre':
nw = 40 if nw is None else nw
freqs, wts = _get_scaled_legendre_roots(nw)
elif alg == 'minimax':
nw = 30 if nw is None else nw
E_max = 0
E_min = 1e10
for k in range(self.nkpts):
for s in range(2):
E_max = max(E_max, mo_energy[s, k, -1] - mo_energy[s, k, 0])
E_min = min(E_min, mo_energy[s, k, self.nocc[s]] - mo_energy[s, k, self.nocc[s] - 1])
E_range = E_max / E_min
freqs, wts = get_frequency_minimax_grid(nw, E_range)
else:
raise NotImplementedError('Grids algorithm not implemented!')
return freqs, wts
[docs]
def get_q_mesh(self, mo_energy, mo_coeff):
"""Get q-mesh for finite size correction.
Equation 39-42 in doi.org/10.1021/acs.jctc.0c00704
Parameters
----------
mo_energy : double 3d array
orbital energy
mo_coeff : double 4d array
coefficient from AO to MO
Returns
-------
qij : double 1d array
q-mesh grids
q_abs : double 1d array
absolute positions of q-mesh grids
nq_pts : init
number of q-mesh grids
"""
# Set up q mesh for q->0 finite size correction
nmoa, nmob = self.nmo
nocca, noccb = self.nocc
nkpts = self.nkpts
if not self.fc_grid:
q_pts = np.array([1e-3, 0, 0]).reshape(1, 3)
else:
Nq = 4
q_pts = np.zeros((Nq**3 - 1, 3))
for i in range(Nq):
for j in range(Nq):
for k in range(Nq):
if i == 0 and j == 0 and k == 0:
continue
else:
q_pts[i * Nq**2 + j * Nq + k - 1, 0] = k * 5e-4
q_pts[i * Nq**2 + j * Nq + k - 1, 1] = j * 5e-4
q_pts[i * Nq**2 + j * Nq + k - 1, 2] = i * 5e-4
nq_pts = len(q_pts)
q_abs = self.mol.get_abs_kpts(q_pts)
# Get qij = 1/sqrt(Omega) * < psi_{ik} | e^{iqr} | psi_{ak-q} > at q: (nkpts, nocc, nvir)
qij_a = np.zeros((nq_pts, nkpts, nocca, nmoa - nocca), dtype=np.complex128)
qij_b = np.zeros((nq_pts, nkpts, noccb, nmob - noccb), dtype=np.complex128)
for k in range(nq_pts):
qij_tmp = get_qij(rpa, q_abs[k], mo_energy, mo_coeff)
qij_a[k] = qij_tmp[0]
qij_b[k] = qij_tmp[1]
return qij_a, qij_b, q_abs, nq_pts
if __name__ == '__main__':
from pyscf.pbc import gto, dft, scf
from pyscf.pbc.lib import chkfile
import os
cell = gto.Cell()
cell.build(
unit='B',
a=[[0.0, 6.74027466, 6.74027466], [6.74027466, 0.0, 6.74027466], [6.74027466, 6.74027466, 0.0]],
atom="""H 0 0 0
H 1.68506866 1.68506866 1.68506866
H 3.37013733 3.37013733 3.37013733""",
basis='gth-dzvp',
pseudo='gth-pade',
verbose=5,
charge=0,
spin=1,
)
cell.spin = cell.spin * 3
kpts = cell.make_kpts([3, 1, 1], scaled_center=[0, 0, 0])
gdf = df.RSDF(cell, kpts)
gdf_fname = 'h3_ints_311.h5'
gdf._cderi_to_save = gdf_fname
if not os.path.isfile(gdf_fname):
gdf.build()
chkfname = 'h_311.chk'
if os.path.isfile(chkfname):
kmf = scf.KUHF(cell, kpts, exxdiv='ewald').rs_density_fit()
kmf.with_df = gdf
kmf.with_df._cderi = gdf_fname
data = chkfile.load(chkfname, 'scf')
kmf.__dict__.update(data)
else:
kmf = scf.KUHF(cell, kpts, exxdiv='ewald').rs_density_fit()
kmf.with_df = gdf
kmf.with_df._cderi = gdf_fname
kmf.conv_tol = 1e-12
kmf.chkfile = chkfname
kmf.kernel()
rpa = KURPA(kmf)
rpa.fc = False
rpa.kernel()
if rank == 0:
print(rpa.e_tot, rpa.e_corr)
assert abs(rpa.e_corr - -0.04288352903004621) < 1e-6
assert abs(rpa.e_tot - -1.584806462873674) < 1e-6
rpa = KURPA(kmf)
rpa.fc = False
rpa.outcore = True
rpa.segsize = 3
rpa.kernel()
if rank == 0:
print(rpa.e_tot, rpa.e_corr)
assert abs(rpa.e_corr - -0.04288352903004621) < 1e-6
assert abs(rpa.e_tot - -1.584806462873674) < 1e-6
# with finite size corrections
rpa = KURPA(kmf)
rpa.fc = True
rpa.kernel()
if rank == 0:
print(rpa.e_tot, rpa.e_corr)
assert abs(rpa.e_corr - -0.04295466718074476) < 1e-6
# Test on Na (metallic)
cell = gto.Cell()
cell.build(
unit='angstrom',
a="""
-2.11250000000000 2.11250000000000 2.11250000000000
2.11250000000000 -2.11250000000000 2.11250000000000
2.11250000000000 2.11250000000000 -2.11250000000000
""",
atom="""
Na 0.00000 0.00000 0.00000
""",
dimension=3,
max_memory=126000,
verbose=5,
pseudo='gth-pade',
basis='gth-dzvp-molopt-sr',
precision=1e-10,
)
kpts = cell.make_kpts([2, 2, 1], scaled_center=[0, 0, 0])
gdf = df.RSDF(cell, kpts)
gdf_fname = 'gdf_ints_221.h5'
gdf._cderi_to_save = gdf_fname
if not os.path.isfile(gdf_fname):
gdf.build()
chkfname = 'na_221.chk'
if os.path.isfile(chkfname):
kmf = dft.KUKS(cell, kpts).rs_density_fit()
kmf = scf.addons.smearing_(kmf, sigma=5e-3, method='fermi')
kmf.xc = 'lda'
kmf.with_df = gdf
kmf.with_df._cderi = gdf_fname
data = chkfile.load(chkfname, 'scf')
kmf.__dict__.update(data)
else:
kmf = dft.KUKS(cell, kpts).rs_density_fit()
kmf = scf.addons.smearing_(kmf, sigma=5e-3, method='fermi')
kmf.xc = 'lda'
kmf.with_df = gdf
kmf.with_df._cderi = gdf_fname
kmf.conv_tol = 1e-12
kmf.chkfile = chkfname
kmf.kernel()
rpa = KURPA(kmf)
rpa.kernel()
assert abs(rpa.e_corr - -0.04031792880477294) < 1e-6
assert abs(rpa.e_tot - -47.60693294927457) < 1e-6
rpa = KURPA(kmf)
rpa.outcore = True
rpa.segsize = 3
rpa.kernel()
assert abs(rpa.e_corr - -0.04031792880477294) < 1e-6
assert abs(rpa.e_tot - -47.60693294927457) < 1e-6
print('passed tests!')