def _add_vvvv_full(mycc, t1T, t2T, eris, out=None, with_ovvv=False):
    '''Ht2 = numpy.einsum('ijcd,acdb->ijab', t2, vvvv)
    without using symmetry t2[ijab] = t2[jiba] in t2 or Ht2
    '''
    time0 = time.clock(), time.time()
    log = logger.Logger(mycc.stdout, mycc.verbose)

    nvir_seg, nvir, nocc = t2T.shape[:3]
    vloc0, vloc1 = _task_location(nvir, rank)
    nocc2 = nocc*(nocc+1)//2
    if t1T is None:
        tau = lib.pack_tril(t2T.reshape(nvir_seg*nvir,nocc2))
    else:
        tau = t2T + numpy.einsum('ai,bj->abij', t1T[vloc0:vloc1], t1T)
        tau = lib.pack_tril(tau.reshape(nvir_seg*nvir,nocc2))
    tau = tau.reshape(nvir_seg,nvir,nocc2)

    if mycc.direct:   # AO-direct CCSD
        if with_ovvv:
            raise NotImplementedError
        mo = getattr(eris, 'mo_coeff', None)
        if mo is None:  # If eris does not have the attribute mo_coeff
            mo = _mo_without_core(mycc, mycc.mo_coeff)

        ao_loc = mycc.mol.ao_loc_nr()
        nao, nmo = mo.shape
        ntasks = mpi.pool.size
        task_sh_locs = lib.misc._balanced_partition(ao_loc, ntasks)
        ao_loc0 = ao_loc[task_sh_locs[rank  ]]
        ao_loc1 = ao_loc[task_sh_locs[rank+1]]

        orbv = mo[:,nocc:]
        tau = lib.einsum('abij,pb->apij', tau, orbv)
        tau_priv = numpy.zeros((ao_loc1-ao_loc0,nao,nocc,nocc))
        for task_id, tau in _rotate_tensor_block(tau):
            loc0, loc1 = _task_location(nvir, task_id)
            tau_priv += lib.einsum('pa,abij->pbij', orbv[ao_loc0:ao_loc1,loc0:loc1], tau)
        tau = None
        time1 = log.timer_debug1('vvvv-tau mo2ao', *time0)

        buf = _contract_vvvv_t2(mycc, None, tau_priv, task_sh_locs, None, log)
        buf = buf.reshape(tau_priv.shape)
        tau_priv = None
        time1 = log.timer_debug1('vvvv-tau contraction', *time1)

        buf = lib.einsum('apij,pb->abij', buf, orbv)
        Ht2 = numpy.ndarray(t2T.shape, buffer=out)
        Ht2[:] = 0
        for task_id, buf in _rotate_tensor_block(buf):
            ao_loc0 = ao_loc[task_sh_locs[task_id  ]]
            ao_loc1 = ao_loc[task_sh_locs[task_id+1]]
            Ht2 += lib.einsum('pa,pbij->abij', orbv[ao_loc0:ao_loc1,vloc0:vloc1], buf)

        time1 = log.timer_debug1('vvvv-tau ao2mo', *time1)
    else:
        raise NotImplementedError
    return Ht2.reshape(t2T.shape)

def amplitudes_to_vector(t1, t2, out=None):
    t2T = t2.transpose(2,3,0,1)
    nvir_seg, nvir, nocc = t2T.shape[:3]
    if rank == 0:
        t1T = t1.T
        nov = nocc * nvir
        nocc2 = nocc*(nocc+1)//2
        size = nov + nvir_seg*nvir*nocc2
        vector = numpy.ndarray(size, t1.dtype, buffer=out)
        vector[:nov] = t1T.ravel()
        lib.pack_tril(t2T.reshape(nvir_seg*nvir,nocc,nocc), out=vector[nov:])
    else:
        vector = lib.pack_tril(t2T.reshape(nvir_seg*nvir,nocc,nocc))
    return vector

    def __init__(self, myci, mo_coeff, method='incore'):
        mol = myci.mol
        mf = myci._scf
        nocc = myci.nocc
        nmo = myci.nmo
        nvir = nmo - nocc
        if mo_coeff is None:
            self.mo_coeff = mo_coeff = myci.mo_coeff
        if (method == 'incore' and mf._eri is not None):
            eri = ao2mo.kernel(mf._eri, mo_coeff, verbose=myci.verbose)
        else:
            eri = ao2mo.kernel(mol, mo_coeff, verbose=myci.verbose)
        eri = ao2mo.restore(1, eri, nmo)
        eri = eri.reshape(nmo,nmo,nmo,nmo)

        self.oooo = eri[:nocc,:nocc,:nocc,:nocc]
        self.vvoo = eri[nocc:,nocc:,:nocc,:nocc]
        self.vooo = eri[nocc:,:nocc,:nocc,:nocc]
        self.voov = eri[nocc:,:nocc,:nocc,nocc:]
        self.vovv = lib.pack_tril(eri[nocc:,:nocc,nocc:,nocc:].reshape(-1,nvir,nvir))
        self.vvvv = ao2mo.restore(4, eri[nocc:,nocc:,nocc:,nocc:].copy(), nvir)

        dm = mf.make_rdm1()
        vhf = mf.get_veff(mol, dm)
        h1 = mf.get_hcore(mol)
        self.fock = reduce(numpy.dot, (mo_coeff.T, h1 + vhf, mo_coeff))

def ecp_int(cell, kpts=None):
    if rank == 0:
        comm.bcast(cell.dumps())
    else:
        cell = pgto.loads(comm.bcast(None))

    if kpts is None:
        kpts_lst = numpy.zeros((1,3))
    else:
        kpts_lst = numpy.reshape(kpts, (-1,3))

    ecpcell = gto.Mole()
    ecpcell._atm = cell._atm
    # append a fictitious s function to mimic the auxiliary index in pbc.incore.
    # ptr2last_env_idx to force PBCnr3c_fill_* function to copy the entire "env"
    ptr2last_env_idx = len(cell._env) - 1
    ecpbas = numpy.vstack([[0, 0, 1, 1, 0, ptr2last_env_idx, 0, 0],
                           cell._ecpbas]).astype(numpy.int32)
    ecpcell._bas = ecpbas
    ecpcell._env = cell._env
    # In pbc.incore _ecpbas is appended to two sets of cell._bas and the
    # fictitious s function.
    cell._env[AS_ECPBAS_OFFSET] = cell.nbas * 2 + 1
    cell._env[AS_NECPBAS] = len(cell._ecpbas)

    kptij_lst = numpy.hstack((kpts_lst,kpts_lst)).reshape(-1,2,3)
    nkpts = len(kpts_lst)
    if abs(kpts_lst).sum() < 1e-9:  # gamma_point
        dtype = numpy.double
    else:
        dtype = numpy.complex128
    ao_loc = cell.ao_loc_nr()
    nao = ao_loc[-1]
    mat = numpy.zeros((nkpts,nao,nao), dtype=dtype)

    intor = cell._add_suffix('ECPscalar')
    int3c = incore.wrap_int3c(cell, ecpcell, intor, kptij_lst=kptij_lst)

    # shls_slice of auxiliary index (0,1) corresponds to the fictitious s function
    tasks = [(i, i+1, j, j+1, 0, 1) # shls_slice
             for i in range(cell.nbas) for j in range(i+1)]
    for shls_slice in mpi.work_stealing_partition(tasks):
        i0 = ao_loc[shls_slice[0]]
        i1 = ao_loc[shls_slice[1]]
        j0 = ao_loc[shls_slice[2]]
        j1 = ao_loc[shls_slice[3]]
        buf = numpy.empty((nkpts,i1-i0,j1-j0), dtype=dtype)
        mat[:,i0:i1,j0:j1] = int3c(shls_slice, buf)

    buf = mpi.reduce(mat)
    if rank == 0:
        mat = []
        for k, kpt in enumerate(kpts_lst):
            v = lib.unpack_tril(lib.pack_tril(buf[k]), lib.HERMITIAN)
            if abs(kpt).sum() < 1e-9:  # gamma_point:
                v = v.real
            mat.append(v)
        if kpts is None or numpy.shape(kpts) == (3,):
            mat = mat[0]
        return mat

def cosmo_fock_o1(cosmo, dm):
    mol = cosmo.mol
    nao = dm.shape[0]
    # phi
    cosmo.loadsegs()
    coords = cosmo.cosurf[:cosmo.nps*3].reshape(-1,3)
    fakemol = _make_fakemol(coords)
    j3c = df.incore.aux_e2(mol, fakemol, intor='cint3c2e_sph', aosym='s2ij')
    tril_dm = lib.pack_tril(dm) * 2
    diagidx = numpy.arange(nao)
    diagidx = diagidx*(diagidx+1)//2 + diagidx
    tril_dm[diagidx] *= .5
    cosmo.phi = -numpy.einsum('x,xk->k', tril_dm, j3c)
    for ia in range(mol.natm):
        cosmo.phi += mol.atom_charge(ia)/lib.norm(mol.atom_coord(ia)-coords, axis=1)
    cosmo.savesegs()
    # qk
    cosmo.charges()
    # vpot
    cosmo.loadsegs()
#X    fakemol = _make_fakemol(cosmo.cosurf[:cosmo.nps*3].reshape(-1,3))
#X    j3c = df.incore.aux_e2(mol, fakemol, intor='cint3c2e_sph', aosym='s2ij')
    fock = lib.unpack_tril(numpy.einsum('xk,k->x', j3c, -cosmo.qcos[:cosmo.nps]))
    fepsi = cosmo.fepsi() 
    fock = fepsi*fock
    return fock

def part_eri_hermi(eri, norb, nimp):
    eri1 = ao2mo.restore(4, eri, norb)
    for i in range(eri1.shape[0]):
        tmp = lib.unpack_tril(eri1[i])
        tmp[nimp:] = 0
        eri1[i] = lib.pack_tril(tmp + tmp.T)
    eri1 = lib.transpose_sum(eri1, inplace=True)
    return ao2mo.restore(8, eri1, norb) * 0.25

def _make_Lpq(mydf, mol, auxmol):
    atm, bas, env, ao_loc = incore._env_and_aoloc('cint3c1e_sph', mol, auxmol)
    nao = ao_loc[mol.nbas]
    naux = ao_loc[-1] - nao
    nao_pair = nao * (nao+1) // 2

    if mydf.metric.upper() == 'S':
        intor = 'cint3c1e_sph'
        s_aux = auxmol.intor_symmetric('cint1e_ovlp_sph')
    elif mydf.metric.upper() == 'T':
        intor = 'cint3c1e_p2_sph'
        s_aux = auxmol.intor_symmetric('cint1e_kin_sph') * 2
    else:  # metric.upper() == 'J'
        intor = 'cint3c2e_sph'
        s_aux = incore.fill_2c2e(mol, auxmol)
    cintopt = gto.moleintor.make_cintopt(atm, bas, env, intor)

    if mydf.charge_constraint:
        ovlp = lib.pack_tril(mol.intor_symmetric('cint1e_ovlp_sph'))

        aux_loc = auxmol.ao_loc_nr()
        s_index = numpy.hstack([range(aux_loc[i],aux_loc[i+1])
                                for i,l in enumerate(auxmol._bas[:,ANG_OF]) if l == 0])
        a = numpy.zeros((naux+1,naux+1))
        a[:naux,:naux] = s_aux
        a[naux,s_index] = a[s_index,naux] = 1
        try:
            cd = scipy.linalg.cho_factor(a)
            def solve(Lpq):
                return scipy.linalg.cho_solve(cd, Lpq)
        except scipy.linalg.LinAlgError:
            def solve(Lpq):
                return scipy.linalg.solve(a, Lpq)
    else:
        cd = scipy.linalg.cho_factor(s_aux)
        def solve(Lpq):
            return scipy.linalg.cho_solve(cd, Lpq, overwrite_b=True)

    def get_Lpq(shls_slice, col0, col1, buf):
        # Be cautious here, _ri.nr_auxe2 assumes buf in F-order
        Lpq = _ri.nr_auxe2(intor, atm, bas, env, shls_slice, ao_loc,
                           's2ij', 1, cintopt, buf).T
        if mydf.charge_constraint:
            Lpq = numpy.ndarray(shape=(naux+1,col1-col0), buffer=buf)
            Lpq[naux,:] = ovlp[col0:col1]
            Lpq1 = solve(Lpq)
            assert(Lpq1.flags.f_contiguous)
            lib.transpose(Lpq1.T, out=Lpq)
            return Lpq[:naux]
        else:
            return solve(Lpq)
    return get_Lpq

    def test_unpack(self):
        a = numpy.random.random((400,400))
        a = a+a*.5j
        for i in range(400):
            a[i,i] = a[i,i].real
        b = a-a.T.conj()
        b = numpy.array((b,b))
        x = lib.hermi_triu(b[0].T, hermi=2, inplace=0)
        self.assertAlmostEqual(abs(b[0].T-x).max(), 0, 12)

        x = lib.hermi_triu(b[1], hermi=2, inplace=0)
        self.assertAlmostEqual(abs(b[1]-x).max(), 0, 12)
        self.assertAlmostEqual(abs(x - lib.unpack_tril(lib.pack_tril(x), 2)).max(), 0, 12)

        x = lib.hermi_triu(a, hermi=1, inplace=0)
        self.assertAlmostEqual(abs(x-x.T.conj()).max(), 0, 12)

        xs = numpy.asarray((x,x,x))
        self.assertAlmostEqual(abs(xs - lib.unpack_tril(lib.pack_tril(xs))).max(), 0, 12)

        numpy.random.seed(1)
        a = numpy.random.random((5050,20))
        self.assertAlmostEqual(lib.finger(lib.unpack_tril(a, axis=0)), -103.03970592075423, 10)

def _int_nuc_vloc(mydf, nuccell, kpts, intor='int3c2e_sph', aosym='s2', comp=1):
    '''Vnuc - Vloc'''
    cell = mydf.cell
    nkpts = len(kpts)

# Use the 3c2e code with steep s gaussians to mimic nuclear density
    fakenuc = aft._fake_nuc(cell)
    fakenuc._atm, fakenuc._bas, fakenuc._env = \
            gto.conc_env(nuccell._atm, nuccell._bas, nuccell._env,
                         fakenuc._atm, fakenuc._bas, fakenuc._env)

    kptij_lst = numpy.hstack((kpts,kpts)).reshape(-1,2,3)
    ishs = mpi.work_balanced_partition(numpy.arange(cell.nbas),
                                       costs=numpy.arange(1, cell.nbas+1))
    if len(ishs) > 0:
        ish0, ish1 = ishs[0], ishs[-1]+1
        buf = incore.aux_e2(cell, fakenuc, intor, aosym='s2', kptij_lst=kptij_lst,
                            shls_slice=(ish0,ish1,0,cell.nbas,0,fakenuc.nbas))
    else:
        buf = numpy.zeros(0)

    charge = cell.atom_charges()
    charge = numpy.append(charge, -charge)  # (charge-of-nuccell, charge-of-fakenuc)
    nao = cell.nao_nr()
    nchg = len(charge)
    nao_pair = nao*(nao+1)//2
    buf = buf.reshape(nkpts,-1,nchg)
# scaled by 1./mpi.pool.size because nuc is mpi.reduced in get_nuc function
    buf = numpy.einsum('kxz,z->kx', buf, 1./mpi.pool.size*charge)
    mat = numpy.empty((nkpts,nao_pair), dtype=numpy.complex128)
    for k in range(nkpts):
        mat[k] = mpi.allgather(buf[k])

    if (rank == 0 and
        cell.dimension == 3 and intor in ('int3c2e', 'int3c2e_sph',
                                          'int3c2e_cart')):
        assert(comp == 1)
        charges = cell.atom_charges()

        nucbar = sum([z/nuccell.bas_exp(i)[0] for i,z in enumerate(charges)])
        nucbar *= numpy.pi/cell.vol

        ovlp = cell.pbc_intor('int1e_ovlp_sph', 1, lib.HERMITIAN, kpts)
        for k in range(nkpts):
            if aosym == 's1':
                mat[k] += nucbar * ovlp[k].reshape(nao_pair)
            else:
                mat[k] += nucbar * lib.pack_tril(ovlp[k])

    return mat

def contract_1e(f1e, fcivec, norb, nelec, link_index=None):
    fcivec = numpy.asarray(fcivec, order='C')
    link_index = _unpack(norb, nelec, link_index)
    na, nlink = link_index.shape[:2]
    assert(fcivec.size == na**2)
    ci1 = numpy.empty_like(fcivec)
    f1e_tril = lib.pack_tril(f1e)
    libfci.FCIcontract_1e_spin0(f1e_tril.ctypes.data_as(ctypes.c_void_p),
                                fcivec.ctypes.data_as(ctypes.c_void_p),
                                ci1.ctypes.data_as(ctypes.c_void_p),
                                ctypes.c_int(norb), ctypes.c_int(na),
                                ctypes.c_int(nlink),
                                link_index.ctypes.data_as(ctypes.c_void_p))
# no *.5 because FCIcontract_2e_spin0 only compute half of the contraction
    return lib.transpose_sum(ci1, inplace=True).reshape(fcivec.shape)

def _int_nuc_vloc(mydf, nuccell, kpts, intor='int3c2e', aosym='s2', comp=1):
    '''Vnuc - Vloc'''
    cell = mydf.cell
    nkpts = len(kpts)

# Use the 3c2e code with steep s gaussians to mimic nuclear density
    fakenuc = _fake_nuc(cell)
    fakenuc._atm, fakenuc._bas, fakenuc._env = \
            gto.conc_env(nuccell._atm, nuccell._bas, nuccell._env,
                         fakenuc._atm, fakenuc._bas, fakenuc._env)

    kptij_lst = numpy.hstack((kpts,kpts)).reshape(-1,2,3)
    buf = incore.aux_e2(cell, fakenuc, intor, aosym=aosym, comp=comp,
                        kptij_lst=kptij_lst)

    charge = cell.atom_charges()
    charge = numpy.append(charge, -charge)  # (charge-of-nuccell, charge-of-fakenuc)
    nao = cell.nao_nr()
    nchg = len(charge)
    if aosym == 's1':
        nao_pair = nao**2
    else:
        nao_pair = nao*(nao+1)//2
    if comp == 1:
        buf = buf.reshape(nkpts,nao_pair,nchg)
        mat = numpy.einsum('kxz,z->kx', buf, charge)
    else:
        buf = buf.reshape(nkpts,comp,nao_pair,nchg)
        mat = numpy.einsum('kcxz,z->kcx', buf, charge)

    # vbar is the interaction between the background charge
    # and the compensating function.  0D, 1D, 2D do not have vbar.
    if cell.dimension == 3 and intor in ('int3c2e', 'int3c2e_sph',
                                         'int3c2e_cart'):
        assert(comp == 1)
        charge = -cell.atom_charges()

        nucbar = sum([z/nuccell.bas_exp(i)[0] for i,z in enumerate(charge)])
        nucbar *= numpy.pi/cell.vol

        ovlp = cell.pbc_intor('int1e_ovlp', 1, lib.HERMITIAN, kpts)
        for k in range(nkpts):
            if aosym == 's1':
                mat[k] -= nucbar * ovlp[k].reshape(nao_pair)
            else:
                mat[k] -= nucbar * lib.pack_tril(ovlp[k])
    return mat

def incore(eri, dm, hermi=0):
    assert(not numpy.iscomplexobj(eri))
    eri = numpy.ascontiguousarray(eri)
    dm = numpy.ascontiguousarray(dm)
    nao = dm.shape[0]
    vj = numpy.empty((nao,nao))
    vk = numpy.empty((nao,nao))
    npair = nao*(nao+1)//2
    if eri.ndim == 2 and npair*npair == eri.size: # 4-fold symmetry eri
        fdrv = getattr(libcvhf, 'CVHFnrs4_incore_drv')
        # 'ijkl,kl->ij'
        fvj = _fpointer('CVHFics4_kl_s2ij')
        # 'ijkl,il->jk'
        fvk = _fpointer('CVHFics4_il_s1jk')
        # or
        ## 'ijkl,ij->kl'
        #fvj = _fpointer('CVHFics4_ij_s2kl')
        ## 'ijkl,jk->il'
        #fvk = _fpointer('CVHFics4_jk_s1il')

        tridm = dm
    elif eri.ndim == 1 and npair*(npair+1)//2 == eri.size: # 8-fold symmetry eri
        fdrv = getattr(libcvhf, 'CVHFnrs8_incore_drv')
        fvj = _fpointer('CVHFics8_tridm_vj')
        if hermi == 1:
            fvk = _fpointer('CVHFics8_jk_s2il')
        else:
            fvk = _fpointer('CVHFics8_jk_s1il')
        tridm = lib.pack_tril(lib.transpose_sum(dm))
        i = numpy.arange(nao)
        tridm[i*(i+1)//2+i] *= .5
    else:
        raise RuntimeError('Array shape not consistent: DM %s, eri %s'
                           % (dm.shape, eri.shape))
    fdrv(eri.ctypes.data_as(ctypes.c_void_p),
         tridm.ctypes.data_as(ctypes.c_void_p),
         vj.ctypes.data_as(ctypes.c_void_p),
         dm.ctypes.data_as(ctypes.c_void_p),
         vk.ctypes.data_as(ctypes.c_void_p),
         ctypes.c_int(nao), fvj, fvk)
    if hermi != 0:
        vj = lib.hermi_triu(vj, hermi)
        vk = lib.hermi_triu(vk, hermi)
    else:
        vj = lib.hermi_triu(vj, 1)
    return vj, vk

def make_phi(pcmobj, dm, r_vdw, ui):
    mol = pcmobj.mol
    natm = mol.natm
    coords_1sph, weights_1sph = make_grids_one_sphere(pcmobj.lebedev_order)
    ngrid_1sph = coords_1sph.shape[0]

    if not (isinstance(dm, numpy.ndarray) and dm.ndim == 2):
        dm = dm[0] + dm[1]
    tril_dm = lib.pack_tril(dm+dm.T)
    nao = dm.shape[0]
    diagidx = numpy.arange(nao)
    diagidx = diagidx*(diagidx+1)//2 + diagidx
    tril_dm[diagidx] *= .5

    atom_coords = mol.atom_coords()
    atom_charges = mol.atom_charges()

    extern_point_idx = ui > 0
    cav_coords = (atom_coords.reshape(natm,1,3)
                  + numpy.einsum('r,gx->rgx', r_vdw, coords_1sph))

    v_phi = numpy.empty((natm,ngrid_1sph))
    for ia in range(natm):
# Note (-) sign is not applied to atom_charges, because (-) is explicitly
# included in rhs and L matrix
        d_rs = atom_coords.reshape(-1,1,3) - cav_coords[ia]
        v_phi[ia] = numpy.einsum('z,zp->p', atom_charges, 1./lib.norm(d_rs,axis=2))

    max_memory = pcmobj.max_memory - lib.current_memory()[0]
    blksize = int(max(max_memory*1e6/8/nao**2, 400))

    cav_coords = cav_coords[extern_point_idx]
    v_phi_e = numpy.empty(cav_coords.shape[0])
    int3c2e = mol._add_suffix('int3c2e')
    for i0, i1 in lib.prange(0, cav_coords.shape[0], blksize):
        fakemol = gto.fakemol_for_charges(cav_coords[i0:i1])
        v_nj = df.incore.aux_e2(mol, fakemol, intor=int3c2e, aosym='s2ij')
        v_phi_e[i0:i1] = numpy.einsum('x,xk->k', tril_dm, v_nj)
    v_phi[extern_point_idx] -= v_phi_e

    ylm_1sph = numpy.vstack(sph.real_sph_vec(coords_1sph, pcmobj.lmax, True))
    phi = -numpy.einsum('n,xn,jn,jn->jx', weights_1sph, ylm_1sph, ui, v_phi)
    return phi

def cosmo_occ_o1(cosmo, dm):
    mol = cosmo.mol
    nao = dm.shape[0]
    #cosmo.check()
    cosmo.occ0()
    cosmo.loadsegs()
    #cosmo.check()
    ioff = 3*cosmo.nps
    coords = cosmo.cosurf[ioff:ioff+cosmo.npspher*3].reshape(-1,3)
    fakemol = _make_fakemol(coords)
    j3c = df.incore.aux_e2(mol, fakemol, intor='cint3c2e_sph', aosym='s2ij')
    tril_dm = lib.pack_tril(dm) * 2
    diagidx = numpy.arange(nao)
    diagidx = diagidx*(diagidx+1)//2 + diagidx
    tril_dm[diagidx] *= .5
    cosmo.phio = -numpy.einsum('x,xk->k', tril_dm, j3c)
    for ia in range(mol.natm):
        cosmo.phio += mol.atom_charge(ia)/lib.norm(mol.atom_coord(ia)-coords, axis=1)
    cosmo.savesegs()
    return cosmo.occ1()

def _int_nuc_vloc(mydf, nuccell, kpts, intor='cint3c2e_sph'):
    '''Vnuc - Vloc'''
    cell = mydf.cell
    rcut = max(cell.rcut, nuccell.rcut)
    Ls = cell.get_lattice_Ls(rcut=rcut)
    expLk = numpy.asarray(numpy.exp(1j*numpy.dot(Ls, kpts.T)), order='C')
    nkpts = len(kpts)

# Use the 3c2e code with steep s gaussians to mimic nuclear density
    fakenuc = _fake_nuc(cell)
    fakenuc._atm, fakenuc._bas, fakenuc._env = \
            gto.conc_env(nuccell._atm, nuccell._bas, nuccell._env,
                         fakenuc._atm, fakenuc._bas, fakenuc._env)

    nao = cell.nao_nr()
    buf = [numpy.zeros((nao,nao,fakenuc.natm), order='F', dtype=numpy.complex128)
           for k in range(nkpts)]
    ints = incore._wrap_int3c(cell, fakenuc, intor, 1, Ls, buf)
    atm, bas, env = ints._envs[:3]
    c_shls_slice = (ctypes.c_int*6)(0, cell.nbas, cell.nbas, cell.nbas*2,
                                    cell.nbas*2, cell.nbas*2+fakenuc.natm)

    xyz = numpy.asarray(cell.atom_coords(), order='C')
    ptr_coordL = atm[:cell.natm,gto.PTR_COORD]
    ptr_coordL = numpy.vstack((ptr_coordL,ptr_coordL+1,ptr_coordL+2)).T.copy('C')
    for l, L1 in enumerate(Ls):
        env[ptr_coordL] = xyz + L1
        exp_Lk = numpy.einsum('k,ik->ik', expLk[l].conj(), expLk[:l+1])
        exp_Lk = numpy.asarray(exp_Lk, order='C')
        exp_Lk[l] = .5
        ints(exp_Lk, c_shls_slice)

    charge = cell.atom_charges()
    charge = numpy.append(charge, -charge)  # (charge-of-nuccell, charge-of-fakenuc)
    for k, kpt in enumerate(kpts):
        v = numpy.einsum('ijz,z->ij', buf[k], charge)
        buf[k] = lib.pack_tril(v + v.T.conj())
    return buf

    def __init__(self, myci, mo_coeff=None, method='incore'):
        cput0 = (time.clock(), time.time())
        moidx = numpy.ones(myci.mo_occ.size, dtype=numpy.bool)
        if isinstance(myci.frozen, (int, numpy.integer)):
            moidx[:myci.frozen] = False
        elif len(myci.frozen) > 0:
            moidx[numpy.asarray(myci.frozen)] = False
        if mo_coeff is None:
            self.mo_coeff = mo_coeff = myci.mo_coeff[:,moidx]
        else:
            self.mo_coeff = mo_coeff = mo_coeff[:,moidx]
        dm = myci._scf.make_rdm1(myci.mo_coeff, myci.mo_occ)
        fockao = myci._scf.get_hcore() + myci._scf.get_veff(myci.mol, dm)
        self.fock = reduce(numpy.dot, (mo_coeff.T, fockao, mo_coeff))

        nocc = myci.nocc
        nmo = myci.nmo
        nvir = nmo - nocc
        mem_incore, mem_outcore, mem_basic = ccsd._mem_usage(nocc, nvir)
        mem_now = lib.current_memory()[0]

        log = logger.Logger(myci.stdout, myci.verbose)
        if (method == 'incore' and myci._scf._eri is not None and
            (mem_incore+mem_now < myci.max_memory) or myci.mol.incore_anyway):
            eri1 = ao2mo.incore.full(myci._scf._eri, mo_coeff)
            #:eri1 = ao2mo.restore(1, eri1, nmo)
            #:self.oooo = eri1[:nocc,:nocc,:nocc,:nocc].copy()
            #:self.ooov = eri1[:nocc,:nocc,:nocc,nocc:].copy()
            #:self.vooo = eri1[nocc:,:nocc,:nocc,:nocc].copy()
            #:self.voov = eri1[nocc:,:nocc,:nocc,nocc:].copy()
            #:self.vvoo = eri1[nocc:,nocc:,:nocc,:nocc].copy()
            #:vovv = eri1[nocc:,:nocc,nocc:,nocc:].copy()
            #:self.vovv = lib.pack_tril(vovv.reshape(-1,nvir,nvir))
            #:self.vvvv = ao2mo.restore(4, eri1[nocc:,nocc:,nocc:,nocc:], nvir)
            nvir_pair = nvir * (nvir+1) // 2
            self.oooo = numpy.empty((nocc,nocc,nocc,nocc))
            self.ooov = numpy.empty((nocc,nocc,nocc,nvir))
            self.vooo = numpy.empty((nvir,nocc,nocc,nocc))
            self.voov = numpy.empty((nvir,nocc,nocc,nvir))
            self.vovv = numpy.empty((nvir,nocc,nvir_pair))
            self.vvvv = numpy.empty((nvir_pair,nvir_pair))
            ij = 0
            outbuf = numpy.empty((nmo,nmo,nmo))
            oovv = numpy.empty((nocc,nocc,nvir,nvir))
            for i in range(nocc):
                buf = lib.unpack_tril(eri1[ij:ij+i+1], out=outbuf[:i+1])
                for j in range(i+1):
                    self.oooo[i,j] = self.oooo[j,i] = buf[j,:nocc,:nocc]
                    self.ooov[i,j] = self.ooov[j,i] = buf[j,:nocc,nocc:]
                    oovv[i,j] = oovv[j,i] = buf[j,nocc:,nocc:]
                ij += i + 1
            self.vvoo = lib.transpose(oovv.reshape(nocc**2,-1)).reshape(nvir,nvir,nocc,nocc)
            oovv = None
            ij1 = 0
            for i in range(nocc,nmo):
                buf = lib.unpack_tril(eri1[ij:ij+i+1], out=outbuf[:i+1])
                self.vooo[i-nocc] = buf[:nocc,:nocc,:nocc]
                self.voov[i-nocc] = buf[:nocc,:nocc,nocc:]
                lib.pack_tril(_cp(buf[:nocc,nocc:,nocc:]), out=self.vovv[i-nocc])
                dij = i - nocc + 1
                lib.pack_tril(_cp(buf[nocc:i+1,nocc:,nocc:]),
                              out=self.vvvv[ij1:ij1+dij])
                ij += i + 1
                ij1 += dij
        else:
            cput1 = time.clock(), time.time()
            _tmpfile1 = tempfile.NamedTemporaryFile(dir=lib.param.TMPDIR)
            _tmpfile2 = tempfile.NamedTemporaryFile(dir=lib.param.TMPDIR)
            self.feri1 = feri1 = h5py.File(_tmpfile1.name)
            def __del__feri1(self):
                feri1.close()
            self.feri1.__del__ = __del__feri1
            orbo = mo_coeff[:,:nocc]
            orbv = mo_coeff[:,nocc:]
            nvpair = nvir * (nvir+1) // 2
            self.oooo = self.feri1.create_dataset('oooo', (nocc,nocc,nocc,nocc), 'f8')
            self.ooov = self.feri1.create_dataset('ooov', (nocc,nocc,nocc,nvir), 'f8')
            self.vvoo = self.feri1.create_dataset('vvoo', (nvir,nvir,nocc,nocc), 'f8')
            self.vooo = self.feri1.create_dataset('vooo', (nvir,nocc,nocc,nocc), 'f8')
            self.voov = self.feri1.create_dataset('voov', (nvir,nocc,nocc,nvir), 'f8')
            self.vovv = self.feri1.create_dataset('vovv', (nvir,nocc,nvpair), 'f8')
            fsort = _ccsd.libcc.CCsd_sort_inplace
            nocc_pair = nocc*(nocc+1)//2
            nvir_pair = nvir*(nvir+1)//2
            def sort_inplace(p0, p1, eri):
                fsort(eri.ctypes.data_as(ctypes.c_void_p),
                      ctypes.c_int(nocc), ctypes.c_int(nvir),
                      ctypes.c_int((p1-p0)*nocc))
                vv = eri[:,:nvir_pair]
                oo = eri[:,nvir_pair:nvir_pair+nocc_pair]
                ov = eri[:,nvir_pair+nocc_pair:].reshape(-1,nocc,nvir)
                return oo, ov, vv
            buf = numpy.empty((nmo,nmo,nmo))
            oovv = numpy.empty((nocc,nocc,nvir,nvir))
            def save_occ_frac(p0, p1, eri):
                oo, ov, vv = sort_inplace(p0, p1, eri)
                self.oooo[p0:p1] = lib.unpack_tril(oo, out=buf).reshape(p1-p0,nocc,nocc,nocc)
                self.ooov[p0:p1] = ov.reshape(p1-p0,nocc,nocc,nvir)
                oovv[p0:p1] = lib.unpack_tril(vv, out=buf).reshape(p1-p0,nocc,nvir,nvir)
            def save_vir_frac(p0, p1, eri):
#.........这里部分代码省略.........

#.........这里部分代码省略.........

    def save1(piece,buf):
        start = piece*chunks_half[1]
        stop = (piece+1)*chunks_half[1]
        if stop > nao_pair:
            stop = nao_pair
        half_eri[:,start:stop] = buf[:,:stop-start]
        return

    def load2(piece):
        start = piece*chunks_half[0]
        stop = (piece+1)*chunks_half[0]
        if stop > nij_pair:
            stop = nij_pair
            if start >= nij_pair:
                start = stop - 1
        return half_eri[start:stop,:]

    def prefetch2(piece):
        start = piece*chunks_half[0]
        stop = (piece+1)*chunks_half[0]
        if stop > nij_pair:
            stop = nij_pair
            if start >= nij_pair:
                start = stop - 1
        buf_prefetch[:stop-start,:] = half_eri[start:stop,:]
        return

    def save2(piece,buf):
        start = piece*chunks_full[0]
        stop = (piece+1)*chunks_full[0]
        if stop > nij_pair:
            stop = nij_pair
        h5d_eri[start:stop,:] = buf[:stop-start,:]
        return

    # transform \mu\nu -> ij
    cput0 = time.clock(), time.time()
    Cimu = mo_coeffs[0].conj().transpose()
    buf_write = numpy.empty((nij_pair,chunks_half[1]))
    buf_out = numpy.empty_like(buf_write)
    wpiece = 0
    with lib.call_in_background(save1) as async_write:
        for lo in range(nao_pair):
            if lo % chunks_half[1] == 0 and lo > 0:
                #save1(wpiece,buf_write)
                buf_out, buf_write = buf_write, buf_out
                async_write(wpiece,buf_out)
                wpiece += 1
            buf = lib.unpack_row(eri,lo)
            uv = lib.unpack_tril(buf)
            uv = Cimu.dot(uv).dot(mo_coeffs[1])
            if ij_red:
                ij = numpy.ravel(uv) # grabs by row
            else:
                ij = lib.pack_tril(uv)
            buf_write[:,lo % chunks_half[1]] = ij
    # final write operation & cleanup
    save1(wpiece,buf_write)
    log.timer('(uv|lo) -> (ij|lo)', *cput0)
    uv = None
    ij = None
    buf = None

    # transform \lambda\sigma -> kl
    cput1 = time.clock(), time.time()
    Cklam = mo_coeffs[2].conj().transpose()
    buf_write = numpy.empty((chunks_full[0],nkl_pair))
    buf_out = numpy.empty_like(buf_write)
    buf_read = numpy.empty((chunks_half[0],nao_pair))
    buf_prefetch = numpy.empty_like(buf_read)
    rpiece = 0
    wpiece = 0
    with lib.call_in_background(save2,prefetch2) as (async_write,prefetch):
        buf_read = load2(rpiece)
        prefetch(rpiece+1)
        for ij in range(nij_pair):
            if ij % chunks_full[0] == 0 and ij > 0:
                #save2(wpiece,buf_write)
                buf_out, buf_write = buf_write, buf_out
                async_write(wpiece,buf_out)
                wpiece += 1
            if ij % chunks_half[0] == 0 and ij > 0:
                #buf_read = load2(rpiece)
                buf_read, buf_prefetch = buf_prefetch, buf_read
                rpiece += 1
                prefetch(rpiece+1)
            lo = lib.unpack_tril(buf_read[ij % chunks_half[0],:])
            lo = Cklam.dot(lo).dot(mo_coeffs[3])
            if kl_red:
                kl = numpy.ravel(lo)
            else:
                kl = lib.pack_tril(lo)
            buf_write[ij % chunks_full[0],:] = kl
    save2(wpiece,buf_write)
    log.timer('(ij|lo) -> (ij|kl)', *cput1)

    if isinstance(erifile, str):
        feri.close()
    return erifile

#.........这里部分代码省略.........

    tau2 = _ccsd.make_tau(t2, t1, t1)
    #:goovv += numpy.einsum('ijkl,klab->ijab', goooo[:,:,j0:j1], tau2)
    lib.dot(goooo.reshape(nocc*nocc,-1),
            tau2.reshape(-1,nvir**2), 1, goovv.reshape(-1,nvir**2), 1)
    tau2 += numpy.einsum('ia,jb->ijab', t1, t1)
    tau2p = tau2.reshape(nocc,nvir,nocc,nvir)
    for i in range(nocc):
        tau2p[i] = tau2[i].transpose(2,0,1)
    tau2, tau2p = tau2p.reshape(nov,-1), None
    #:goovv += numpy.einsum('ibld,jlda->ijab', mOvOv, tau2) * .5
    #:goovv -= numpy.einsum('iald,jldb->ijab', mOVov, tau2) * .5
    tmp = lib.dot(mOvOv.reshape(-1,nov), tau2.T, .5).reshape(nocc,nvir,-1,nvir)
    for i in range(nocc):
        tmp[i] = goovv[i].transpose(1,0,2) + tmp[i].transpose(2,1,0)
    goovv, tmp = tmp, None
    lib.dot(mOVov.reshape(-1,nov), tau2.T, -.5, goovv.reshape(nov,-1), 1)

    #:goovv += numpy.einsum('iald,jlbd->ijab', mOVov*2+mOvOv, t2) * .5
    t2a, tau2 = tau2.reshape(nocc,nvir,nocc,nvir), None
    for i in range(nocc):
        t2a[i] = t2[i].transpose(1,0,2)
    tmp = mOVov*2
    tmp += mOvOv
    lib.dot(tmp.reshape(-1,nov), t2a.reshape(nov,-1), .5, goovv.reshape(nov,-1), 1)
    t2a = tmp = None
    for i in range(nocc):
        goovv[i] = goovv[i] * 2 - goovv[i].transpose(2,1,0)
    dovov = goovv
    goooo = goovv = None

    #:gvovv += numpy.einsum('aick,kb->aibc', pvOVo, t1)
    mOVov = lib.transpose(mOVov.reshape(nov,-1))
    gvovv = lib.dot(mOVov.reshape(nocc,-1).T, t1).reshape(nvir,nocc,nvir,nvir)
    mOVov = None
    tmp = numpy.einsum('ja,jb->ab', l1, t1)
    #:gvovv += numpy.einsum('ab,ic->aibc', tmp, t1)
    #:gvovv += numpy.einsum('ba,ic->aibc', mvv, t1*.5)
    for i in range(nvir):
        gvovv[i] += numpy.einsum('b,ic->icb', tmp[i], t1)
        gvovv[i] += numpy.einsum('b,ic->icb', mvv[:,i]*.5, t1)
        gvovv[i] = gvovv[i].transpose(0,2,1)

    #:gvovv += numpy.einsum('ja,jibc->aibc', l1, t2)
    #:gvovv += numpy.einsum('jibc,ja->aibc', l2, t1)
    #:gvovv -= numpy.einsum('aibk,kc->aibc', pvOvO, t1)
    mOvOv = lib.transpose(mOvOv.reshape(nov,-1))
    lib.dot(mOvOv.reshape(nocc,-1).T, t1, -1, gvovv.reshape(-1,nvir), 1)
    mOvOv = None
    lib.dot(l1.T, t2.reshape(nocc,-1), 1, gvovv.reshape(nvir,-1), 1)
    lib.dot(t1.T, l2.reshape(nocc,-1), 1, gvovv.reshape(nvir,-1), 1)

    tmp = numpy.empty((nocc,nvir,nvir))
    for i in range(nvir):
        #:gvovv*2 - gvovv.transpose(0,1,3,2)
        gvovv[i] = _ccsd.make_021(gvovv[i], gvovv[i], 2, -1, out=tmp)

    #:gvvvv = numpy.einsum('ijab,ijcd->abcd', l2, t2)*.5
    #:jabc = numpy.einsum('ijab,ic->jabc', l2, t1) * .5
    #:gvvvv += numpy.einsum('jabc,jd->abcd', jabc, t1)
    #:gvovv -= numpy.einsum('adbc,id->aibc', gvvvv, t1*2)
    tau = _ccsd.make_tau(t2, t1, t1)
    theta = make_theta(tau)
    tau = None
    l2tmp = lib.pack_tril(l2.reshape(-1,nvir,nvir))
    gtmp = lib.dot(l2tmp.T, theta.reshape(nocc**2,-1), .5).reshape(-1,nvir,nvir)
    l2tmp = theta = None
    nvir_pair = nvir * (nvir+1) //2
    tmp = numpy.empty((nvir,nvir,nvir))
    tmp1 = numpy.empty((nvir,nvir,nvir))
    tmptril = numpy.empty((nvir,nvir_pair))
    diag_idx = numpy.arange(nvir)
    diag_idx = diag_idx*(diag_idx+1)//2 + diag_idx

    dvvvv = numpy.empty((nvir_pair,nvir_pair))
    dovvv = numpy.empty((nocc,nvir,nvir,nvir))
# dvvov = (gvovv*2 - gvovv.transpose(0,1,3,2)).transpose(0,2,1,3)
# dovvv = dvvov.transpose(2,3,0,1)
    p0 = 0
    for i in range(nvir):
        tmp[:i+1] = gtmp[p0:p0+i+1]
        for j in range(i+1, nvir):
            tmp[j] = gtmp[j*(j+1)//2+i].T
        lib.dot(t1, tmp.reshape(nvir,-1), -2, gvovv[i].reshape(nocc,-1), 1)
        dovvv[:,:,i] = gvovv[i].transpose(0,2,1)

        #:gvvvv[i] = (tmp*2-tmp.transpose(0,2,1)).transpose(1,0,2)
        #:gvvvv = .5*(gvvvv+gvvvv.transpose(0,1,3,2))
        #:dvvvv = .5*(gvvvv+gvvvv.transpose(1,0,3,2))
        tmp1[:] = tmp.transpose(1,0,2)
        _ccsd.precontract(tmp1, diag_fac=2, out=tmptril)
        dvvvv[p0:p0+i] += tmptril[:i]
        dvvvv[p0:p0+i] *= .25
        dvvvv[i*(i+1)//2+i] = tmptril[i] * .5
        for j in range(i+1, nvir):
            dvvvv[j*(j+1)//2+i] = tmptril[j]
        p0 += i + 1
    gtmp = tmp = tmp1 = tmptril = gvovv = None
    dvvov = dovvv.transpose(2,3,0,1)
    return (dovov, dvvvv, doooo, doovv, dovvo, dvvov, dovvv, dooov)

    def ft_fuse(job_id, uniq_kptji_id, sh0, sh1):
        kpt = uniq_kpts[uniq_kptji_id]  # kpt = kptj - kpti
        adapted_ji_idx = numpy.where(uniq_inverse == uniq_kptji_id)[0]
        adapted_kptjs = kptjs[adapted_ji_idx]
        nkptj = len(adapted_kptjs)

        Gaux = ft_ao.ft_ao(fused_cell, Gv, None, b, gxyz, Gvbase, kpt).T
        Gaux = fuse(Gaux)
        Gaux *= mydf.weighted_coulG(kpt, False, mesh)
        kLR = lib.transpose(numpy.asarray(Gaux.real, order='C'))
        kLI = lib.transpose(numpy.asarray(Gaux.imag, order='C'))
        j2c = numpy.asarray(fswap['j2c/%d'%uniq_kptji_id])
        j2ctag = j2ctags[uniq_kptji_id]
        naux0 = j2c.shape[0]
        if ('j2c-/%d' % uniq_kptji_id) in fswap:
            j2c_negative = numpy.asarray(fswap['j2c-/%d'%uniq_kptji_id])
        else:
            j2c_negative = None

        if is_zero(kpt):
            aosym = 's2'
        else:
            aosym = 's1'

        if aosym == 's2' and cell.dimension == 3:
            vbar = fuse(mydf.auxbar(fused_cell))
            ovlp = cell.pbc_intor('int1e_ovlp', hermi=1, kpts=adapted_kptjs)
            ovlp = [lib.pack_tril(s) for s in ovlp]

        j3cR = [None] * nkptj
        j3cI = [None] * nkptj
        i0 = ao_loc[sh0]
        i1 = ao_loc[sh1]
        for k, idx in enumerate(adapted_ji_idx):
            key = 'j3c-chunks/%d/%d' % (job_id, idx)
            v = fuse(numpy.asarray(fswap[key]))
            if aosym == 's2' and cell.dimension == 3:
                for i in numpy.where(vbar != 0)[0]:
                    v[i] -= vbar[i] * ovlp[k][i0*(i0+1)//2:i1*(i1+1)//2].ravel()
            j3cR[k] = numpy.asarray(v.real, order='C')
            if v.dtype == numpy.complex128:
                j3cI[k] = numpy.asarray(v.imag, order='C')
            v = None

        ncol = j3cR[0].shape[1]
        Gblksize = max(16, int(max_memory*1e6/16/ncol/(nkptj+1)))  # +1 for pqkRbuf/pqkIbuf
        Gblksize = min(Gblksize, ngrids, 16384)
        pqkRbuf = numpy.empty(ncol*Gblksize)
        pqkIbuf = numpy.empty(ncol*Gblksize)
        buf = numpy.empty(nkptj*ncol*Gblksize, dtype=numpy.complex128)
        log.alldebug2('    blksize (%d,%d)', Gblksize, ncol)

        if aosym == 's2':
            shls_slice = (sh0, sh1, 0, sh1)
        else:
            shls_slice = (sh0, sh1, 0, cell.nbas)
        for p0, p1 in lib.prange(0, ngrids, Gblksize):
            dat = ft_ao._ft_aopair_kpts(cell, Gv[p0:p1], shls_slice, aosym, b,
                                        gxyz[p0:p1], Gvbase, kpt,
                                        adapted_kptjs, out=buf)
            nG = p1 - p0
            for k, ji in enumerate(adapted_ji_idx):
                aoao = dat[k].reshape(nG,ncol)
                pqkR = numpy.ndarray((ncol,nG), buffer=pqkRbuf)
                pqkI = numpy.ndarray((ncol,nG), buffer=pqkIbuf)
                pqkR[:] = aoao.real.T
                pqkI[:] = aoao.imag.T

                lib.dot(kLR[p0:p1].T, pqkR.T, -1, j3cR[k], 1)
                lib.dot(kLI[p0:p1].T, pqkI.T, -1, j3cR[k], 1)
                if not (is_zero(kpt) and gamma_point(adapted_kptjs[k])):
                    lib.dot(kLR[p0:p1].T, pqkI.T, -1, j3cI[k], 1)
                    lib.dot(kLI[p0:p1].T, pqkR.T,  1, j3cI[k], 1)

        for k, idx in enumerate(adapted_ji_idx):
            if is_zero(kpt) and gamma_point(adapted_kptjs[k]):
                v = j3cR[k]
            else:
                v = j3cR[k] + j3cI[k] * 1j
            if j2ctag == 'CD':
                v = scipy.linalg.solve_triangular(j2c, v, lower=True, overwrite_b=True)
                fswap['j3c-chunks/%d/%d'%(job_id,idx)][:naux0] = v
            else:
                fswap['j3c-chunks/%d/%d'%(job_id,idx)][:naux0] = lib.dot(j2c, v)

            # low-dimension systems
            if j2c_negative is not None:
                fswap['j3c-/%d/%d'%(job_id,idx)] = lib.dot(j2c_negative, v)