def h_op(x):
     x = self.unpack_uniq_var(x)
     norb = x.shape[0]
     #:hx = numpy.einsum('qp,xjj,xjq->pj', x, dip, dip)
     #:hx+= numpy.einsum('qp,xqq,xjq->pj', x, dip, dip)
     #:hx+= numpy.einsum('jk,xkk,xkp->pj', x, dip, dip)
     #:hx+= numpy.einsum('jk,xpp,xkp->pj', x, dip, dip)
     #:hx+= numpy.einsum('qj,xjq,xjp->pj', x, dip, dip)
     #:hx+= numpy.einsum('pk,xjp,xkp->pj', x, dip, dip)
     #:hx-= numpy.einsum('qp,xpp,xjq->pj', x, dip, dip) * 2
     #:hx-= numpy.einsum('qp,xjp,xpq->pj', x, dip, dip) * 2
     #:hx+= numpy.einsum('qj,xjp,xjq->pj', x, dip, dip)
     #:hx+= numpy.einsum('pk,xkp,xjp->pj', x, dip, dip)
     #:hx-= numpy.einsum('jk,xjj,xkp->pj', x, dip, dip) * 2
     #:hx-= numpy.einsum('jk,xkj,xjp->pj', x, dip, dip) * 2
     #:return -self.pack_uniq_var(hx)
     #:hx = numpy.einsum('iq,qp->pi', g0, x)
     hx = lib.dot(x.T, g0.T).conj()
     #:hx+= numpy.einsum('qi,xiq,xip->pi', x, dip, dip) * 2
     hx+= numpy.einsum('xip,xi->pi', dip, numpy.einsum('qi,xiq->xi', x, dip)) * 2
     #:hx-= numpy.einsum('qp,xpp,xiq->pi', x, dip, dip) * 2
     hx-= numpy.einsum('xpp,xip->pi', dip,
                       lib.dot(dip.reshape(-1,norb), x).reshape(3,norb,norb)) * 2
     #:hx-= numpy.einsum('qp,xip,xpq->pi', x, dip, dip) * 2
     hx-= numpy.einsum('xip,xp->pi', dip, numpy.einsum('qp,xpq->xp', x, dip)) * 2
     return -self.pack_uniq_var(hx-hx.conj().T)

def general(mydf, mo_coeffs, kpts=None, compact=True):
    if isinstance(mo_coeffs, numpy.ndarray) and mo_coeffs.ndim == 2:
        mo_coeffs = (mo_coeffs,) * 4

    eri = mydf.get_eri(kpts)

####################
# gamma point, the integral is real and with s4 symmetry
    if eri.dtype == numpy.float64:
        return ao2mo.general(eri, mo_coeffs, compact=compact)
    else:
        mokl, klslice = ao2mo.incore._conc_mos(mo_coeffs[2], mo_coeffs[3],
                                               False)[2:]
        if mokl.dtype == numpy.float64:
            mokl = mokl + 0j
        nao = mo_coeffs[0].shape[0]
        nmoi = mo_coeffs[0].shape[1]
        nmoj = mo_coeffs[1].shape[1]
        nmok = mo_coeffs[2].shape[1]
        nmol = mo_coeffs[3].shape[1]
        moi = numpy.asarray(mo_coeffs[0], order='F')
        moj = numpy.asarray(mo_coeffs[1], order='F')
        tao = [0]
        ao_loc = None
        pqkl = _ao2mo.r_e2(eri.reshape(-1,nao**2), mokl, klslice, tao, ao_loc, aosym='s1')
        pqkl = pqkl.reshape(nao,nao,nmok*nmol)
        pjkl = numpy.empty((nao,nmoj,nmok*nmol), dtype=numpy.complex128)
        for i in range(nao):
            lib.dot(moj.T, pqkl[i], 1, pjkl[i], 0)
        pqkl = None
        eri_mo = lib.dot(moi.T.conj(), pjkl.reshape(nao,-1))
        return eri_mo.reshape(nmoi*nmoj,-1)

def get_j_kpts(mf, cell, dm_kpts, kpts, kpt_band=None):
    coords = gen_grid.gen_uniform_grids(cell)
    nkpts = len(kpts)
    ngs = len(coords)
    dm_kpts = np.asarray(dm_kpts)
    nao = dm_kpts.shape[-1]

    ni = numint._KNumInt(kpts)
    aoR_kpts = ni.eval_ao(cell, coords, kpts)
    if kpt_band is not None:
        aoR_kband = numint.eval_ao(cell, coords, kpt_band)

    dms = dm_kpts.reshape(-1,nkpts,nao,nao)
    nset = dms.shape[0]

    vjR = [get_vjR(cell, dms[i], aoR_kpts) for i in range(nset)]
    if kpt_band is not None:
        vj_kpts = [cell.vol/ngs * lib.dot(aoR_kband.T.conj()*vjR[i], aoR_kband)
                   for i in range(nset)]
        if dm_kpts.ndim == 3:  # One set of dm_kpts for KRHF
            vj_kpts = vj_kpts[0]
        return lib.asarray(vj_kpts)
    else:
        vj_kpts = []
        for i in range(nset):
            vj = [cell.vol/ngs * lib.dot(aoR_k.T.conj()*vjR[i], aoR_k)
                  for aoR_k in aoR_kpts]
            vj_kpts.append(lib.asarray(vj))
        return lib.asarray(vj_kpts).reshape(dm_kpts.shape)

def gamma1_intermediates(mycc, t1, t2, l1, l2):
    nocc, nvir = t1.shape
    doo =-numpy.einsum('ja,ia->ij', l1, t1)
    dvv = numpy.einsum('ia,ib->ab', l1, t1)
    dvo = l1.T
    xtv = numpy.einsum('ie,me->im', t1, l1)
    dov = t1 - numpy.einsum('im,ma->ia', xtv, t1)
    #:doo -= numpy.einsum('jkab,ikab->ij', l2, theta)
    #:dvv += numpy.einsum('jica,jicb->ab', l2, theta)
    #:xt1  = numpy.einsum('mnef,inef->mi', l2, make_theta(t2))
    #:xt2  = numpy.einsum('mnaf,mnef->ea', l2, make_theta(t2))
    #:dov += numpy.einsum('imae,me->ia', make_theta(t2), l1)
    #:dov -= numpy.einsum('ma,ie,me->ia', t1, t1, l1)
    #:dov -= numpy.einsum('mi,ma->ia', xt1, t1)
    #:dov -= numpy.einsum('ie,ae->ia', t1, xt2)
    max_memory = mycc.max_memory - lib.current_memory()[0]
    unit = nocc*nvir**2
    blksize = max(ccsd.BLKMIN, int(max_memory*.95e6/8/unit))
    for p0, p1 in prange(0, nocc, blksize):
        theta = make_theta(t2[p0:p1])
        doo[p0:p1] -= lib.dot(theta.reshape(p1-p0,-1), l2.reshape(nocc,-1).T)
        dov[p0:p1] += numpy.einsum('imae,me->ia', theta, l1)
        xt1 = lib.dot(l2.reshape(nocc,-1), theta.reshape(p1-p0,-1).T)
        dov[p0:p1] -= numpy.einsum('mi,ma->ia', xt1, t1)
        xt2 = lib.dot(theta.reshape(-1,nvir).T, l2[p0:p1].reshape(-1,nvir))
        dov -= numpy.einsum('ie,ae->ia', t1, xt2)
        dvv += lib.dot(l2[p0:p1].reshape(-1,nvir).T, theta.reshape(-1,nvir))
    return doo, dov, dvo, dvv

 def get_eri(self):
     nao = self.mol.nao_nr()
     nao_pair = nao * (nao+1) // 2
     ao_eri = numpy.zeros((nao_pair,nao_pair))
     for eri1 in self.loop():
         lib.dot(eri1.T, eri1, 1, ao_eri, 1)
     return ao2mo.restore(8, ao_eri, nao)

 def add_k_(vk, dm, pqk, coulG, buf=None):
     nG = pqk.shape[-1]
     pqk *= numpy.sqrt(coulG)
     ipk = lib.dot(dm, pqk.reshape(nao,-1)).reshape(nao,nao,-1)
     pik = numpy.ndarray((nao,nao,nG), buffer=buf)
     pik[:] = ipk.transpose(1,0,2)
     vk += lib.dot(pqk.reshape(nao,-1), pik.reshape(nao,-1).T)
     return vk

 def cost_function(self, u=None):
     if u is None: u = numpy.eye(self.mo_coeff.shape[1])
     mo_coeff = lib.dot(self.mo_coeff, u)
     dip = dipole_integral(self.mol, mo_coeff)
     r2 = self.mol.intor_symmetric('int1e_r2')
     r2 = numpy.einsum('pi,pi->', mo_coeff, lib.dot(r2, mo_coeff))
     val = r2 - numpy.einsum('xii,xii->', dip, dip) * 2
     return val

def transform_ci_for_orbital_rotation(ci, norb, nelec, u):
    '''Transform CI coefficients to the representation in new one-particle basis.
    Solving CI problem for Hamiltonian h1, h2 defined in old basis,
    CI_old = fci.kernel(h1, h2, ...)
    Given orbital rotation u, the CI problem can be either solved by
    transforming the Hamiltonian, or transforming the coefficients.
    CI_new = fci.kernel(u^T*h1*u, ...) = transform_ci_for_orbital_rotation(CI_old, u)

    Args:
        u : 2D array or a list of 2D array
            the orbital rotation to transform the old one-particle basis to new
            one-particle basis
    '''
    neleca, nelecb = _unpack(nelec)
    strsa = numpy.asarray(cistring.gen_strings4orblist(range(norb), neleca))
    strsb = numpy.asarray(cistring.gen_strings4orblist(range(norb), nelecb))
    one_particle_strs = numpy.asarray([1<<i for i in range(norb)])
    na = len(strsa)
    nb = len(strsb)

    if isinstance(u, numpy.ndarray) and u.ndim == 2:
        ua = ub = u
    else:
        ua, ub = u

    # Unitary transformation array trans_ci is the overlap between two sets of CI basis.
    occ_masks = (strsa[:,None] & one_particle_strs) != 0
    trans_ci_a = numpy.zeros((na,na))
    #for i in range(na): # for old basis
    #    for j in range(na):
    #        uij = u[occ_masks[i]][:,occ_masks[j]]
    #        trans_ci_a[i,j] = numpy.linalg.det(uij)
    occ_idx_all_strs = numpy.where(occ_masks)[1]
    for i in range(na):
        ui = ua[occ_masks[i]].T.copy()
        minors = numpy.take(ui, occ_idx_all_strs, axis=0).reshape(na,neleca,neleca)
        trans_ci_a[i,:] = numpy.linalg.det(minors)

    if neleca == nelecb and numpy.allclose(ua, ub):
        trans_ci_b = trans_ci_a
    else:
        occ_masks = (strsb[:,None] & one_particle_strs) != 0
        trans_ci_b = numpy.zeros((nb,nb))
        #for i in range(nb):
        #    for j in range(nb):
        #        uij = u[occ_masks[i]][:,occ_masks[j]]
        #        trans_ci_b[i,j] = numpy.linalg.det(uij)
        occ_idx_all_strs = numpy.where(occ_masks)[1]
        for i in range(nb):
            ui = ub[occ_masks[i]].T.copy()
            minors = numpy.take(ui, occ_idx_all_strs, axis=0).reshape(nb,nelecb,nelecb)
            trans_ci_b[i,:] = numpy.linalg.det(minors)

    # Transform old basis to new basis for all alpha-electron excitations
    ci = lib.dot(trans_ci_a.T, ci.reshape(na,nb))
    # Transform old basis to new basis for all beta-electron excitations
    ci = lib.dot(ci.reshape(na,nb), trans_ci_b)
    return ci

 def h_op(x):
     x = self.unpack_uniq_var(x)
     norb = x.shape[0]
     hx = lib.dot(x.T, g0.T)
     hx+= numpy.einsum('xip,xi->pi', pop, numpy.einsum('qi,xiq->xi', x, pop)) * 2
     hx-= numpy.einsum('xpp,xip->pi', pop,
                       lib.dot(pop.reshape(-1,norb), x).reshape(-1,norb,norb)) * 2
     hx-= numpy.einsum('xip,xp->pi', pop, numpy.einsum('qp,xpq->xp', x, pop)) * 2
     return -self.pack_uniq_var(hx-hx.T)

def gamma1_intermediates(mycc, t1, t2, l1, l2):
    nocc, nvir = t1.shape
    goo = -numpy.einsum('ja,ia->ij', l1, t1)
    gvv = numpy.einsum('ia,ib->ab', l1, t1)
    #:goo -= numpy.einsum('jkab,ikab->ij', l2, theta)
    #:gvv += numpy.einsum('jica,jicb->ab', l2, theta)
    theta = make_theta(t2)
    goo -= lib.dot(theta.reshape(nocc,-1), l2.reshape(nocc,-1).T)
    gvv += lib.dot(l2.reshape(-1,nvir).T, theta.reshape(-1,nvir))
    return goo, gvv

 def short_range_k(dm, Lpq, j3c):
     Lpq = Lpq.reshape(-1,nao)
     j3c = j3c.reshape(-1,nao)
     iLp = lib.dot(dm, Lpq.T).reshape(-1,nao)
     vk  = lib.dot(j3c.T, iLp)
     if hermi:
         vk = vk + vk.T
     else:
         iLp = lib.dot(dm.T, Lpq.T).reshape(-1,nao)
         vk += lib.dot(iLp.T, j3c)
     return vk

 def shift_grids(r):
     r_frac = lib.dot(r - origin, b.T)
     # Grids on the boundary (r_frac == +/-0.5) of the new cell may lead to
     # unbalanced contributions to the dipole moment. Exclude them from the
     # dipole and quadrupole
     r_frac[r_frac== 0.5] = 0
     r_frac[r_frac==-0.5] = 0
     r_frac[r_frac > 0.5] -= 1
     r_frac[r_frac <-0.5] += 1
     r = lib.dot(r_frac, a)
     return r

 def trans(aoiR, aojR, fac=1):
     if id(aoiR) == id(aojR) and ao2mo.incore.iden_coeffs(mo_coeffs[0], mo_coeffs[1]):
         moiR = mojR = numpy.asarray(lib.dot(mo_coeffs[0].T,aoiR.T), order='C')
     else:
         moiR = numpy.asarray(lib.dot(mo_coeffs[0].T, aoiR.T), order='C')
         mojR = numpy.asarray(lib.dot(mo_coeffs[1].T, aojR.T), order='C')
     mo_pairs_G = numpy.empty((nmoi,nmoj,ngs), dtype=numpy.complex128)
     for i in range(nmoi):
         mo_pairs_G[i] = tools.fft(fac * moiR[i].conj() * mojR, mydf.gs)
     mo_pairs_G = mo_pairs_G.reshape(-1,ngs).T
     return mo_pairs_G

 def trans(aoi, aoj, fac=1):
     if id(aoi) == id(aoj) and iden_coeffs(mo_coeffs[0], mo_coeffs[1]):
         moi = moj = numpy.asarray(lib.dot(mo_coeffs[0].T,aoi.T), order='C')
     else:
         moi = numpy.asarray(lib.dot(mo_coeffs[0].T, aoi.T), order='C')
         moj = numpy.asarray(lib.dot(mo_coeffs[1].T, aoj.T), order='C')
     mo_pairs_G = numpy.empty((nmoi,nmoj,ngrids), dtype=numpy.complex128)
     for i in range(nmoi):
         mo_pairs_G[i] = tools.fft(fac * moi[i].conj() * moj, mydf.mesh)
     mo_pairs_G = mo_pairs_G.reshape(-1,ngrids).T
     return mo_pairs_G

def general(mydf, mo_coeffs, kpts=None,
            compact=getattr(__config__, 'pbc_df_ao2mo_general_compact', True)):
    '''General MO integral transformation'''
    from pyscf.pbc.df.df_ao2mo import warn_pbc2d_eri
    warn_pbc2d_eri(mydf)
    cell = mydf.cell
    kptijkl = _format_kpts(kpts)
    kpti, kptj, kptk, kptl = kptijkl
    if isinstance(mo_coeffs, numpy.ndarray) and mo_coeffs.ndim == 2:
        mo_coeffs = (mo_coeffs,) * 4
    mo_coeffs = [numpy.asarray(mo, order='F') for mo in mo_coeffs]
    if not _iskconserv(cell, kptijkl):
        lib.logger.warn(cell, 'fft_ao2mo: momentum conservation not found in '
                        'the given k-points %s', kptijkl)
        return numpy.zeros([mo.shape[1] for mo in mo_coeffs])

    allreal = not any(numpy.iscomplexobj(mo) for mo in mo_coeffs)
    q = kptj - kpti
    coulG = tools.get_coulG(cell, q, mesh=mydf.mesh)
    coords = cell.gen_uniform_grids(mydf.mesh)
    max_memory = mydf.max_memory - lib.current_memory()[0]

    if gamma_point(kptijkl) and allreal:
        ao = mydf._numint.eval_ao(cell, coords, kpti)[0]
        if ((iden_coeffs(mo_coeffs[0], mo_coeffs[1]) and
             iden_coeffs(mo_coeffs[0], mo_coeffs[2]) and
             iden_coeffs(mo_coeffs[0], mo_coeffs[3]))):
            moiT = mojT = numpy.asarray(lib.dot(mo_coeffs[0].T,ao.T), order='C')
            ao = None
            max_memory = max_memory - moiT.nbytes*1e-6
            eri = _contract_compact(mydf, (moiT,mojT), coulG, max_memory=max_memory)
            if not compact:
                nmo = moiT.shape[0]
                eri = ao2mo.restore(1, eri, nmo).reshape(nmo**2,nmo**2)
        else:
            mos = [numpy.asarray(lib.dot(c.T, ao.T), order='C') for c in mo_coeffs]
            ao = None
            fac = numpy.array(1.)
            max_memory = max_memory - sum([x.nbytes for x in mos])*1e-6
            eri = _contract_plain(mydf, mos, coulG, fac, max_memory=max_memory).real
        return eri

    else:
        aos = mydf._numint.eval_ao(cell, coords, kptijkl)
        mos = [numpy.asarray(lib.dot(c.T, aos[i].T), order='C')
               for i,c in enumerate(mo_coeffs)]
        aos = None
        fac = numpy.exp(-1j * numpy.dot(coords, q))
        max_memory = max_memory - sum([x.nbytes for x in mos])*1e-6
        eri = _contract_plain(mydf, mos, coulG, fac, max_memory=max_memory)
        return eri

    def make_kpt(kpt):  # kpt = kptj - kpti
        # search for all possible ki and kj that has ki-kj+kpt=0
        kk_match = numpy.einsum('ijx->ij', abs(kk_table + kpt)) < 1e-9
        kpti_idx, kptj_idx = numpy.where(kk_todo & kk_match)
        nkptj = len(kptj_idx)
        log.debug1('kpt = %s', kpt)
        log.debug1('kpti_idx = %s', kpti_idx)
        log.debug1('kptj_idx = %s', kptj_idx)
        kk_todo[kpti_idx,kptj_idx] = False
        if swap_2e and abs(kpt).sum() > 1e-9:
            kk_todo[kptj_idx,kpti_idx] = False

        mydf.exxdiv = exxdiv
        vkcoulG = tools.get_coulG(cell, kpt, True, mydf, mydf.gs) / cell.vol
        # <r|-G+k_rs|s> = conj(<s|G-k_rs|r>) = conj(<s|G+k_sr|r>)
        for k, pqkR, pqkI, p0, p1 \
                in mydf.ft_loop(cell, mydf.gs, kpt, kpts[kptj_idx],
                                max_memory=max_memory):
            ki = kpti_idx[k]
            kj = kptj_idx[k]
            coulG = numpy.sqrt(vkcoulG[p0:p1])

# case 1: k_pq = (pi|iq)
            pqkR *= coulG
            pqkI *= coulG
            rsk =(pqkR.reshape(nao,nao,-1).transpose(1,0,2) -
                  pqkI.reshape(nao,nao,-1).transpose(1,0,2)*1j)
            qpk = rsk.conj()
            for i in range(nset):
                qsk = lib.dot(dms[i,kj], rsk.reshape(nao,-1)).reshape(nao,nao,-1)
                #:vk_kpts[i,ki] += numpy.einsum('qpk,qsk->ps', qpk, qsk)
                vk_kpts[i,ki] += lib.dot(qpk.transpose(1,0,2).reshape(nao,-1),
                                         qsk.transpose(1,0,2).reshape(nao,-1).T)
                qsk = None
            rsk = qpk = None

# case 2: k_pq = (iq|pi)
            if swap_2e and abs(kpt).sum() > 1e-9:
                srk = pqkR - pqkI*1j
                pqk = srk.reshape(nao,nao,-1).conj()
                for i in range(nset):
                    prk = lib.dot(dms[i,ki].T, srk.reshape(nao,-1)).reshape(nao,nao,-1)
                    #:vk_kpts[i,kj] += numpy.einsum('prk,pqk->rq', prk, pqk)
                    vk_kpts[i,kj] += lib.dot(prk.transpose(1,0,2).reshape(nao,-1),
                                             pqk.transpose(1,0,2).reshape(nao,-1).T)
                    prk = None
                srk = pqk = None

        pqkR = pqkI = coulG = None
        return None

def gamma1_intermediates(mycc, t1, t2, l1, l2, max_memory=2000):
    nocc, nvir = t1.shape
    goo = -numpy.einsum('ja,ia->ij', l1, t1)
    gvv = numpy.einsum('ia,ib->ab', l1, t1)
    #:goo -= numpy.einsum('jkab,ikab->ij', l2, theta)
    #:gvv += numpy.einsum('jica,jicb->ab', l2, theta)
    max_memory = max_memory - lib.current_memory()[0]
    unit = nocc*nvir**2
    blksize = max(ccsd.BLKMIN, int(max_memory*.95e6/8/unit))
    for p0, p1 in prange(0, nocc, blksize):
        theta = make_theta(t2[p0:p1])
        goo[p0:p1] -= lib.dot(theta.reshape(p1-p0,-1), l2.reshape(nocc,-1).T)
        gvv += lib.dot(l2[p0:p1].reshape(-1,nvir).T, theta.reshape(-1,nvir))
    return goo, gvv

def _ifftn_blas(g, mesh):
    Gx = np.fft.fftfreq(mesh[0])
    Gy = np.fft.fftfreq(mesh[1])
    Gz = np.fft.fftfreq(mesh[2])
    expRGx = np.exp(np.einsum('x,k->xk', 2j*np.pi*np.arange(mesh[0]), Gx))
    expRGy = np.exp(np.einsum('x,k->xk', 2j*np.pi*np.arange(mesh[1]), Gy))
    expRGz = np.exp(np.einsum('x,k->xk', 2j*np.pi*np.arange(mesh[2]), Gz))
    out = np.empty(g.shape, dtype=np.complex128)
    buf = np.empty(mesh, dtype=np.complex128)
    for i, gi in enumerate(g):
        buf[:] = gi.reshape(mesh)
        f = lib.dot(buf.reshape(mesh[0],-1).T, expRGx, 1./mesh[0], c=out[i].reshape(-1,mesh[0]))
        f = lib.dot(f.reshape(mesh[1],-1).T, expRGy, 1./mesh[1], c=buf.reshape(-1,mesh[1]))
        f = lib.dot(f.reshape(mesh[2],-1).T, expRGz, 1./mesh[2], c=out[i].reshape(-1,mesh[2]))
    return out.reshape(-1, *mesh)

def _fftn_blas(f, mesh):
    Gx = np.fft.fftfreq(mesh[0])
    Gy = np.fft.fftfreq(mesh[1])
    Gz = np.fft.fftfreq(mesh[2])
    expRGx = np.exp(np.einsum('x,k->xk', -2j*np.pi*np.arange(mesh[0]), Gx))
    expRGy = np.exp(np.einsum('x,k->xk', -2j*np.pi*np.arange(mesh[1]), Gy))
    expRGz = np.exp(np.einsum('x,k->xk', -2j*np.pi*np.arange(mesh[2]), Gz))
    out = np.empty(f.shape, dtype=np.complex128)
    buf = np.empty(mesh, dtype=np.complex128)
    for i, fi in enumerate(f):
        buf[:] = fi.reshape(mesh)
        g = lib.dot(buf.reshape(mesh[0],-1).T, expRGx, c=out[i].reshape(-1,mesh[0]))
        g = lib.dot(g.reshape(mesh[1],-1).T, expRGy, c=buf.reshape(-1,mesh[1]))
        g = lib.dot(g.reshape(mesh[2],-1).T, expRGz, c=out[i].reshape(-1,mesh[2]))
    return out.reshape(-1, *mesh)

def get_eri(mydf, kpts=None, compact=False):
    cell = mydf.cell
    if kpts is None:
        kptijkl = numpy.zeros((4,3))
    elif numpy.shape(kpts) == (3,):
        kptijkl = numpy.vstack([kpts]*4)
    else:
        kptijkl = numpy.reshape(kpts, (4,3))

    kpti, kptj, kptk, kptl = kptijkl
    nao = cell.nao_nr()
    nao_pair = nao * (nao+1) // 2
    coulG = tools.get_coulG(cell, kptj-kpti, gs=mydf.gs)
    ngs = len(coulG)

####################
# gamma point, the integral is real and with s4 symmetry
    if abs(kptijkl).sum() < 1e-9:
        ao_pairs_G = get_ao_pairs_G(mydf, kptijkl[:2], compact)
        ao_pairs_G *= numpy.sqrt(coulG).reshape(-1,1)
        aoijR = ao_pairs_G.real.copy()
        aoijI = ao_pairs_G.imag.copy()
        ao_pairs_G = None
        eri = lib.dot(aoijR.T, aoijR, cell.vol/ngs**2)
        eri = lib.dot(aoijI.T, aoijI, cell.vol/ngs**2, eri, 1)
        return eri

####################
# (kpt) i == j == k == l != 0
# (kpt) i == l && j == k && i != j && j != k  =>
#
# complex integrals, N^4 elements
    elif (abs(kpti-kptl).sum() < 1e-9) and (abs(kptj-kptk).sum() < 1e-9):
        ao_pairs_G = get_ao_pairs_G(mydf, kptijkl[:2], False)
        ao_pairs_G *= numpy.sqrt(coulG).reshape(-1,1)
        ao_pairs_invG = ao_pairs_G.T.reshape(nao,nao,-1).transpose(1,0,2).conj()
        ao_pairs_invG = ao_pairs_invG.reshape(-1,ngs)
        return lib.dot(ao_pairs_G.T, ao_pairs_invG.T, cell.vol/ngs**2)

####################
# aosym = s1, complex integrals
#
    else:
        ao_pairs_G = get_ao_pairs_G(mydf, kptijkl[:2], False)
# ao_pairs_invG = rho_rs(-G+k_rs) = conj(rho_sr(G+k_sr)).swap(r,s)
        ao_pairs_invG = get_ao_pairs_G(mydf, -kptijkl[2:], False).conj()
        ao_pairs_G *= coulG.reshape(-1,1)
        return lib.dot(ao_pairs_G.T, ao_pairs_invG, cell.vol/ngs**2)