def exdiagH(h1e, g2e, norb, nelec, writefile=True):
    '''
        exactly diagonalize the hamiltonian.
    '''
    h2e = direct_spin1.absorb_h1e(h1e, g2e, norb, nelec, .5)
    if isinstance(nelec, (int, np.integer)):
        nelecb = nelec//2
        neleca = nelec - nelecb
    else:
        neleca, nelecb = nelec

    naa = cistring.num_strings(norb, neleca)
    nbb = cistring.num_strings(norb, nelecb)
    ndim = naa*nbb
    eyebas = np.eye(ndim, ndim)
    def hop(c):
        hc = direct_spin1.contract_2e(h2e, c, norb, nelec)
        return hc.reshape(-1)
    Hmat = []
    for i in range(ndim):
        hc = hop(eyebas[i].copy())
        Hmat.append(hc)

    Hmat = np.asarray(Hmat)
#    Hmat = Hmat.T.copy()
    ew, ev = nl.eigh(Hmat.T)
    if writefile:
        np.savetxt("cards/eignE.dat", ew, fmt="%10.10f")
        np.savetxt("cards/eignV.dat", ev, fmt="%10.10f")
    return ew, ev

def ft_rdm1s(h1e, g2e, norb, nelec, T, m=50, nsamp=40, Tmin=10e-4):
    '''rdm of spin a and b at temperature T
    '''
    if T < Tmin:
       e, c = kernel(h1e, g2e, norb, nelec)
       rdma, rdmb = direct_spin1.make_rdm1s(c, norb, nelec)
       return rdma, rdmb

    h2e = direct_spin1.absorb_h1e(h1e, g2e, norb, nelec, .5)
    if isinstance(nelec, (int, numpy.integer)):
        nelecb = nelec//2
        neleca = nelec - nelecb
    else:
        neleca, nelecb = nelec
    na = cistring.num_strings(norb, neleca)
    nb = cistring.num_strings(norb, nelecb)
    def vecgen(n1=na, n2=nb):
        ci0 = numpy.random.randn(n1, n2)
#        ci0[0, 0] = 1.
        return ci0.reshape(-1)
    def hop(c):
        hc = direct_spin1.contract_2e(h2e, c, norb, nelec)
        return hc.reshape(-1)
    def qud(v1, v2):
        dma, dmb = direct_spin1.trans_rdm1s(v1, v2, norb, nelec)
        return dma, dmb

#    rdma, rdmb = flan.ht_rdm1s(qud, hop, vecgen, T, norb, m, nsamp)
    rdma, rdmb = flan.ftlan_rdm1s(qud, hop, vecgen, T, norb, m, nsamp)
    return rdma, rdmb

def kernel(h1e, eri, norb, nelec, ecore=0, verbose=logger.NOTE):
    if isinstance(nelec, (int, numpy.integer)):
        nelecb = nelec//2
        neleca = nelec - nelecb
    else:
        neleca, nelecb = nelec
    h2e = direct_spin1.absorb_h1e(h1e, eri, norb, nelec, .5)
    namax = cistring.num_strings(norb, neleca)
    nbmax = cistring.num_strings(norb, nelecb)

    myci = SelectedCI()

    strsa = [int('1'*neleca, 2)]
    strsb = [int('1'*nelecb, 2)]
    ci_strs = (strsa, strsb)
    ci0 = numpy.ones((1,1))
    ci0, ci_strs = enlarge_space(myci, (ci0, ci_strs), h2e, norb, nelec)

    def all_linkstr_index(ci_strs):
        cd_indexa = cre_des_linkstr(ci_strs[0], norb, neleca)
        dd_indexa = des_des_linkstr(ci_strs[0], norb, neleca)
        cd_indexb = cre_des_linkstr(ci_strs[1], norb, nelecb)
        dd_indexb = des_des_linkstr(ci_strs[1], norb, nelecb)
        return cd_indexa, dd_indexa, cd_indexb, dd_indexb

    def hop(c):
        hc = contract_2e(h2e, (c, ci_strs), norb, nelec, link_index)
        return hc.reshape(-1)
    precond = lambda x, e, *args: x/(hdiag-e+1e-4)

    e_last = 0
    tol = 1e-2
    conv = False
    for icycle in range(norb):
        tol = max(tol*1e-2, myci.float_tol)
        link_index = all_linkstr_index(ci_strs)
        hdiag = make_hdiag(h1e, eri, ci_strs, norb, nelec)
        e, ci0 = lib.davidson(hop, ci0.reshape(-1), precond, tol=tol,
                              verbose=verbose)
        print('icycle %d  ci.shape %s  E = %.15g' %
              (icycle, (len(ci_strs[0]), len(ci_strs[1])), e))
        if ci0.shape == (namax,nbmax) or abs(e-e_last) < myci.float_tol*10:
            conv = True
            break
        ci1, ci_strs = enlarge_space(myci, (ci0, ci_strs), h2e, norb, nelec)
        if ci1.size < ci0.size*1.02:
            conv = True
            break
        e_last = e
        ci0 = ci1

    link_index = all_linkstr_index(ci_strs)
    hdiag = make_hdiag(h1e, eri, ci_strs, norb, nelec)
    e, ci0 = lib.davidson(hop, ci0.reshape(-1), precond, tol=myci.conv_tol,
                          verbose=verbose)

    na = len(ci_strs[0])
    nb = len(ci_strs[1])
    return e+ecore, (ci0.reshape(na,nb), ci_strs)

def kernel_fixed_space(myci, h1e, eri, norb, nelec, ci_strs, ci0=None,
                       tol=None, lindep=None, max_cycle=None, max_space=None,
                       nroots=None, davidson_only=None,
                       max_memory=None, verbose=None, ecore=0, **kwargs):
    if verbose is None:
        log = logger.Logger(myci.stdout, myci.verbose)
    elif isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(myci.stdout, verbose)
    if tol is None: tol = myci.conv_tol
    if lindep is None: lindep = myci.lindep
    if max_cycle is None: max_cycle = myci.max_cycle
    if max_space is None: max_space = myci.max_space
    if max_memory is None: max_memory = myci.max_memory
    if nroots is None: nroots = myci.nroots
    if myci.verbose >= logger.WARN:
        myci.check_sanity()

    nelec = direct_spin1._unpack_nelec(nelec, myci.spin)
    ci0, nelec, ci_strs = _unpack(ci0, nelec, ci_strs)
    na = len(ci_strs[0])
    nb = len(ci_strs[1])
    h2e = direct_spin1.absorb_h1e(h1e, eri, norb, nelec, .5)
    h2e = ao2mo.restore(1, h2e, norb)

    link_index = _all_linkstr_index(ci_strs, norb, nelec)
    hdiag = myci.make_hdiag(h1e, eri, ci_strs, norb, nelec)

    if isinstance(ci0, _SCIvector):
        if ci0.size == na*nb:
            ci0 = [ci0.ravel()]
        else:
            ci0 = [x.ravel() for x in ci0]
    else:
        ci0 = myci.get_init_guess(ci_strs, norb, nelec, nroots, hdiag)

    def hop(c):
        hc = myci.contract_2e(h2e, _as_SCIvector(c, ci_strs), norb, nelec, link_index)
        return hc.reshape(-1)
    precond = lambda x, e, *args: x/(hdiag-e+1e-4)

    #e, c = lib.davidson(hop, ci0, precond, tol=myci.conv_tol)
    e, c = myci.eig(hop, ci0, precond, tol=tol, lindep=lindep,
                    max_cycle=max_cycle, max_space=max_space, nroots=nroots,
                    max_memory=max_memory, verbose=log, **kwargs)
    if nroots > 1:
        return e+ecore, [_as_SCIvector(ci.reshape(na,nb),ci_strs) for ci in c]
    else:
        return e+ecore, _as_SCIvector(c.reshape(na,nb), ci_strs)

def kernel(h1e, g2e, norb, nelec):

    h2e = direct_spin1.absorb_h1e(h1e, g2e, norb, nelec, .5)
    if isinstance(nelec, (int, numpy.integer)):
        neleca = nelec//2
    else:
        neleca = nelec[0]

    na = cistring.num_strings(norb, neleca)
    ci0 = numpy.zeros((na,na))
    ci0[0,0] = 1

    def hop(c):
        hc = direct_spin1.contract_2e(h2e, c, norb, nelec)
        return hc.reshape(-1)
    hdiag = direct_spin1.make_hdiag(h1e, g2e, norb, nelec)
    precond = lambda x, e, *args: x/(hdiag-e+1e-4)
    e, c = pyscf.lib.davidson(hop, ci0.reshape(-1), precond)
    return e, c

def kernel_ft_smpl(h1e, g2e, norb, nelec, T, vecgen = 0, m=50, nsmpl = 250, nblk = 10, Tmin=10e-4, nrotation = 200):
    if T < Tmin:
        e, c = kernel(h1e, g2e, norb, nelec)
        return e
    disp = numpy.exp(T) * 0.1 # displacement
    h2e = direct_spin1.absorb_h1e(h1e, g2e, norb, nelec, .5)
    if isinstance(nelec, (int, numpy.integer)):
        nelecb = nelec//2
        neleca = nelec - nelecb
    else:
        neleca, nelecb = nelec

    na = cistring.num_strings(norb, neleca)
    nb = cistring.num_strings(norb, nelecb)
    ci0 = numpy.random.randn(na, nb)
    def hop(c):
        hc = direct_spin1.contract_2e(h2e, c, norb, nelec)
        return hc.reshape(-1)

    E, dev, ar = ftsmpl(hop, ci0, T, flan.ftlan_E1c, nsamp = nsmpl, dr = disp, genci=vecgen, nblock = nblk, nrot = nrotation)
    # ar is the acceptance ratio
    return E, dev, ar

def kernel_ft(h1e, g2e, norb, nelec, T, m=50, nsamp=100, Tmin=10e-4):
    '''E at temperature T
    '''
    if T < Tmin:
        e, c = kernel(h1e, g2e, norb, nelec)
        return e
    h2e = direct_spin1.absorb_h1e(h1e, g2e, norb, nelec, .5)
    if isinstance(nelec, (int, numpy.integer)):
        nelecb = nelec//2
        neleca = nelec - nelecb
    else:
        neleca, nelecb = nelec

    na = cistring.num_strings(norb, neleca)
    nb = cistring.num_strings(norb, nelecb)
    def vecgen(n1=na, n2=nb):
        ci0 = numpy.random.randn(n1, n2)
        return ci0.reshape(-1)
    def hop(c):
        hc = direct_spin1.contract_2e(h2e, c, norb, nelec)
        return hc.reshape(-1)
    E = flan.ftlan_E(hop, vecgen, T, m, nsamp)
    return E

def energy(h1e, eri, fcivec, norb, nelec, link_index=None):
    h2e = direct_spin1.absorb_h1e(h1e, eri, norb, nelec, .5)
    ci1 = contract_2e(h2e, fcivec, norb, nelec, link_index)
    return numpy.dot(fcivec.ravel(), ci1.ravel())

 def absorb_h1e(self, h1e, eri, norb, nelec, fac=1):
     nelec = _unpack_nelec(nelec, self.spin)
     return direct_spin1.absorb_h1e(h1e, eri, norb, nelec, fac)

def energy(h1e, eri, fcivec, norb, nelec, link_index=None, orbsym=None, wfnsym=0):
    h2e = direct_spin1.absorb_h1e(h1e, eri, norb, nelec) * .5
    ci1 = contract_2e(h2e, fcivec, norb, nelec, link_index, orbsym, wfnsym)
    return numpy.dot(fcivec.ravel(), ci1.ravel())

def energy(h1e, eri, fcivec, norb, nelec, link_index=None):
    from pyscf.fci import direct_spin1
    h2e = direct_spin1.absorb_h1e(h1e, eri, norb, nelec, .5)
    ci1 = direct_spin1.contract_2e(h2e, fcivec, norb, nelec, link_index)
    return numpy.dot(fcivec.reshape(-1), ci1.reshape(-1))

    jk = jk.ravel()
    jk_sorted = abs(jk).argsort()[::-1]
    ci1 = [as_SCIvector(numpy.ones(1), hf_str)]

    myci = SelectedCI()
    myci.select_cutoff = .001
    myci.ci_coeff_cutoff = .001

    ci2 = enlarge_space(myci, ci1, h1, eri, jk, eri_sorted, jk_sorted, norb, nelec)
    print(len(ci2[0]))

    ci2 = enlarge_space(myci, ci1, h1, eri, jk, eri_sorted, jk_sorted, norb, nelec)
    numpy.random.seed(1)
    ci3 = numpy.random.random(ci2[0].size)
    ci3 *= 1./numpy.linalg.norm(ci3)
    ci3 = [ci3]
    ci3 = enlarge_space(myci, ci2, h1, eri, jk, eri_sorted, jk_sorted, norb, nelec)

    efci = direct_spin1.kernel(h1, eri, norb, nelec, verbose=5)[0]

    ci4 = contract_2e_ctypes((h1, eri), ci3[0], norb, nelec)

    fci3 = to_fci(ci3, norb, nelec)
    h2e = direct_spin1.absorb_h1e(h1, eri, norb, nelec, .5)
    fci4 = direct_spin1.contract_2e(h2e, fci3, norb, nelec)
    fci4 = from_fci(fci4, ci3[0]._strs, norb, nelec)
    print(abs(ci4-fci4).sum())

    e = myci.kernel(h1, eri, norb, nelec, verbose=5)[0]
    print(e, efci)

 def absorb_h1e(self, h1e, eri, norb, nelec, fac=1):
     return direct_spin1.absorb_h1e(h1e, eri, norb, nelec, fac)

def kernel_float_space(myci, h1e, eri, norb, nelec, ci0=None,
                       tol=None, lindep=None, max_cycle=None, max_space=None,
                       nroots=None, davidson_only=None,
                       max_memory=None, verbose=None, ecore=0, **kwargs):
    if verbose is None:
        log = logger.Logger(myci.stdout, myci.verbose)
    elif isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(myci.stdout, verbose)
    if tol is None: tol = myci.conv_tol
    if lindep is None: lindep = myci.lindep
    if max_cycle is None: max_cycle = myci.max_cycle
    if max_space is None: max_space = myci.max_space
    if max_memory is None: max_memory = myci.max_memory
    if nroots is None: nroots = myci.nroots
    if myci.verbose >= logger.WARN:
        myci.check_sanity()

    nelec = direct_spin1._unpack_nelec(nelec, myci.spin)
    h2e = direct_spin1.absorb_h1e(h1e, eri, norb, nelec, .5)
    h2e = ao2mo.restore(1, h2e, norb)

# TODO: initial guess from CISD
    if isinstance(ci0, _SCIvector):
        if ci0.size == len(ci0._strs[0])*len(ci0._strs[1]):
            ci0 = [ci0.ravel()]
        else:
            ci0 = [x.ravel() for x in ci0]
    else:
        if isinstance(nelec, (int, numpy.integer)):
            nelecb = nelec//2
            neleca = nelec - nelecb
            nelec = neleca, nelecb
        ci_strs = (numpy.asarray([int('1'*nelec[0], 2)]),
                   numpy.asarray([int('1'*nelec[1], 2)]))
        ci0 = _as_SCIvector(numpy.ones((1,1)), ci_strs)
        ci0 = myci.enlarge_space(ci0, h2e, norb, nelec)
        if ci0.size < nroots:
            log.warn('''
  Selected-CI space generated from HF ground state (by double exciting) is not enough for excited states.
  HOMO->LUMO excitations are included in the initial guess.
  NOTE: This may introduce excited states of different symmetry.\n''')
            corea = '1' * (nelec[0]-1)
            coreb = '1' * (nelec[1]-1)
            ci_strs = (numpy.asarray([int('1'+corea, 2), int('10'+corea, 2)]),
                       numpy.asarray([int('1'+coreb, 2), int('10'+coreb, 2)]))
            ci0 = _as_SCIvector(numpy.ones((2,2)), ci_strs)
            ci0 = myci.enlarge_space(ci0, h2e, norb, nelec)
        if ci0.size < nroots:
            raise RuntimeError('Not enough selected-CI space for %d states' % nroots)
        ci_strs = ci0._strs
        hdiag = myci.make_hdiag(h1e, eri, ci_strs, norb, nelec)
        ci0 = myci.get_init_guess(ci_strs, norb, nelec, nroots, hdiag)

    def hop(c):
        hc = myci.contract_2e(h2e, _as_SCIvector(c, ci_strs), norb, nelec, link_index)
        return hc.ravel()
    precond = lambda x, e, *args: x/(hdiag-e+myci.level_shift)

    namax = cistring.num_strings(norb, nelec[0])
    nbmax = cistring.num_strings(norb, nelec[1])
    e_last = 0
    float_tol = 3e-4
    conv = False
    for icycle in range(norb):
        ci_strs = ci0[0]._strs
        float_tol = max(float_tol*.3, tol*1e2)
        log.debug('cycle %d  ci.shape %s  float_tol %g',
                  icycle, (len(ci_strs[0]), len(ci_strs[1])), float_tol)

        ci0 = [c.ravel() for c in ci0]
        link_index = _all_linkstr_index(ci_strs, norb, nelec)
        hdiag = myci.make_hdiag(h1e, eri, ci_strs, norb, nelec)
        #e, ci0 = lib.davidson(hop, ci0.reshape(-1), precond, tol=float_tol)
        e, ci0 = myci.eig(hop, ci0, precond, tol=float_tol, lindep=lindep,
                          max_cycle=max_cycle, max_space=max_space, nroots=nroots,
                          max_memory=max_memory, verbose=log, **kwargs)
        if nroots > 1:
            ci0 = [_as_SCIvector(c, ci_strs) for c in ci0]
            de, e_last = min(e)-e_last, min(e)
            log.info('cycle %d  E = %s  dE = %.8g', icycle, e+ecore, de)
        else:
            ci0 = [_as_SCIvector(ci0, ci_strs)]
            de, e_last = e-e_last, e
            log.info('cycle %d  E = %.15g  dE = %.8g', icycle, e+ecore, de)

        if ci0[0].shape == (namax,nbmax) or abs(de) < tol*1e3:
            conv = True
            break

        last_ci0_size = float(len(ci_strs[0])), float(len(ci_strs[1]))
        ci0 = myci.enlarge_space(ci0, h2e, norb, nelec)
        na = len(ci0[0]._strs[0])
        nb = len(ci0[0]._strs[1])
        if ((.99 < na/last_ci0_size[0] < 1.01) and
            (.99 < nb/last_ci0_size[1] < 1.01)):
            conv = True
            break

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