def METHOD_NAME(self):
    rho0 = kmf.get_rho()
    kpts = libkpts.make_kpts(cell, kpts0, space_group_symmetry=True, time_reversal_symmetry=True)
    dms_ibz = kmf.make_rdm1()[kpts.ibz2bz]
    nao = dms_ibz.shape[-1]
    rho = 0.
    for k in range(kpts.nkpts_ibz):
        rho_k = khf.get_rho(kmf, dms_ibz[k].reshape((-1,nao,nao)), kpts=kpts.kpts_ibz[k].reshape((-1,3)))
        rho += kpts.symmetrize_density(rho_k, k, cell.mesh)
    rho *= 1.0 / kpts.nkpts
    self.assertAlmostEqual(abs(rho - rho0).max(), 0, 8)