def METHOD_NAME(mf, cell=None, dm=None, hermi=0, kpt=None, kpts_band=None,
           with_j=True, with_k=True, **kwargs):
    if cell is None: cell = mf.cell
    if dm is None: dm = mf.make_rdm1()
    if kpt is None: kpt = mf.kpt
    dm = np.asarray(dm)
    nso = dm.shape[-1]
    nao = nso // 2
    dms = dm.reshape(-1,nso,nso)
    if kpts_band is None:  # kpts_band is set to kpt
        nband = 1
        jk_shape = dm.shape
    elif getattr(kpts_band, 'ndim', None) == 1:  # single kpt_band
        nband = 1
        jk_shape = dm.shape
    elif dm.ndim == 2:
        nband = len(kpts_band)
        jk_shape = (nband, nso, nso)
    else:
        nband = len(kpts_band)
        jk_shape = (dm.shape[0], nband, nso, nso)
    dmaa = dms[:,:nao,:nao]
    dmab = dms[:,nao:,:nao]
    dmbb = dms[:,nao:,nao:]
    if with_k:
        if hermi:
            dms = np.stack((dmaa, dmbb, dmab))
        else:
            dmba = dms[:,nao:,:nao]
            dms = np.stack((dmaa, dmbb, dmab, dmba))
        # Note the off-diagonal block breaks the hermitian
        _hermi = 0
    else:
        dms = np.stack((dmaa, dmbb))
        _hermi = 1
    nblocks, n_dm = dms.shape[:2]
    dms = dms.reshape(nblocks*n_dm, nao, nao)
    if mf.rsjk:
        logger.warn(mf, 'RSJK does not support GHF')
    j1, k1 = mf.with_df.METHOD_NAME(dms, _hermi, kpt, kpts_band, with_j, with_k,
                               exxdiv=mf.exxdiv)
    vj = vk = None
    if with_j:
        j1 = j1.reshape(nblocks,n_dm*nband,nao,nao)
        vj = np.zeros((n_dm*nband,nso,nso), j1.dtype)
        vj[:,:nao,:nao] = vj[:,nao:,nao:] = j1[0] + j1[1]
        vj = vj.reshape(jk_shape)
    if with_k:
        k1 = k1.reshape(nblocks,n_dm*nband,nao,nao)
        vk = np.zeros((n_dm*nband,nso,nso), k1.dtype)
        vk[:,:nao,:nao] = k1[0]
        vk[:,nao:,nao:] = k1[1]
        vk[:,:nao,nao:] = k1[2]
        if hermi:
            vk[:,nao:,:nao] = k1[2].conj().transpose(0,2,1)
        else:
            vk[:,nao:,:nao] = k1[3]
        vk = vk.reshape(jk_shape)
    return vj, vk