def get_grad(self, mo_coeff, mo_occ, fock=None):
        mol = self.mol
        if not mol.symmetry:
            return uhf.UHF.get_grad(self, mo_coeff, mo_occ, fock)

        if fock is None:
            dm1 = self.make_rdm1(mo_coeff, mo_occ)
            fock = self.get_hcore(mol) + self.get_veff(self.mol, dm1)
        ovlp_ao = self.get_ovlp()
        orbsyma = symm.label_orb_symm(self, mol.irrep_id, mol.symm_orb,
                                      mo_coeff[0], ovlp_ao, False)
        orbsymb = symm.label_orb_symm(self, mol.irrep_id, mol.symm_orb,
                                      mo_coeff[1], ovlp_ao, False)
        orbsyma = numpy.asarray(orbsyma)
        orbsymb = numpy.asarray(orbsymb)

        occidxa = mo_occ[0] > 0
        occidxb = mo_occ[1] > 0
        viridxa = ~occidxa
        viridxb = ~occidxb
        ga = reduce(numpy.dot, (mo_coeff[0][:,viridxa].T.conj(), fock[0],
                                mo_coeff[0][:,occidxa]))
        ga[orbsyma[viridxa].reshape(-1,1)!=orbsyma[occidxa]] = 0
        gb = reduce(numpy.dot, (mo_coeff[1][:,viridxb].T.conj(), fock[1],
                                mo_coeff[1][:,occidxb]))
        gb[orbsymb[viridxb].reshape(-1,1)!=orbsymb[occidxb]] = 0
        return numpy.hstack((ga.ravel(), gb.ravel()))

def get_irrep_nelec(mol, mo_coeff, mo_occ, s=None):
    '''Alpha/beta electron numbers for each irreducible representation.

    Args:
        mol : an instance of :class:`Mole`
            To provide irrep_id, and spin-adapted basis
        mo_occ : a list of 1D ndarray
            Regular occupancy, without grouping for irreps
        mo_coeff : a list of 2D ndarray
            Regular orbital coefficients, without grouping for irreps

    Returns:
        irrep_nelec : dict
            The number of alpha/beta electrons for each irrep {'ir_name':(int,int), ...}.

    Examples:

    >>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', symmetry=True, charge=1, spin=1, verbose=0)
    >>> mf = scf.UHF(mol)
    >>> mf.scf()
    -75.623975516256721
    >>> scf.uhf_symm.get_irrep_nelec(mol, mf.mo_coeff, mf.mo_occ)
    {'A1': (3, 3), 'A2': (0, 0), 'B1': (1, 1), 'B2': (1, 0)}
    '''
    orbsyma = symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb,
                                  mo_coeff[0], s)
    orbsymb = symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb,
                                  mo_coeff[1], s)
    orbsyma = numpy.array(orbsyma)
    orbsymb = numpy.array(orbsymb)
    irrep_nelec = dict([(mol.irrep_name[k], (int(sum(mo_occ[0][orbsyma==ir])),
                                             int(sum(mo_occ[1][orbsymb==ir]))))
                        for k, ir in enumerate(mol.irrep_id)])
    return irrep_nelec

 def gen_g_hop(self, mo_coeff, mo_occ, fock_ao=None, h1e=None):
     mol = self._scf.mol
     if mol.symmetry:
         ovlp_ao = self._scf.get_ovlp()
         self._orbsym =(symm.label_orb_symm(mol, mol.irrep_id,
                                 mol.symm_orb, mo_coeff[0], s=ovlp_ao),
                        symm.label_orb_symm(mol, mol.irrep_id,
                                 mol.symm_orb, mo_coeff[1], s=ovlp_ao))
     return gen_g_hop_uhf(self, mo_coeff, mo_occ, fock_ao, h1e)

    def kernel(self, mo_coeff=None, ci0=None, macro=None, micro=None, callback=None, _kern=None):
        if mo_coeff is None:
            mo_coeff = self.mo_coeff
        else:
            self.mo_coeff = mo_coeff
        if macro is None:
            macro = self.max_cycle_macro
        if micro is None:
            micro = self.max_cycle_micro
        if callback is None:
            callback = self.callback
        if _kern is None:
            _kern = mc1step.kernel

        if self.verbose > logger.QUIET:
            pyscf.gto.mole.check_sanity(self, self._keys, self.stdout)

        self.mol.check_sanity(self)
        self.dump_flags()

        # irrep_name = self.mol.irrep_name
        irrep_name = self.mol.irrep_id
        try:
            self.orbsym = symm.label_orb_symm(self.mol, irrep_name, self.mol.symm_orb, mo_coeff, s=self._scf.get_ovlp())
        except ValueError:
            logger.warn(self, "mc1step_symm symmetrizes input orbitals")
            s = self._scf.get_ovlp()
            mo_coeff = symm.symmetrize_orb(self.mol, mo_coeff, s=s)
            diag = numpy.einsum("ki,ki->i", mo_coeff, numpy.dot(s, mo_coeff))
            mo_coeff = numpy.einsum("ki,i->ki", mo_coeff, 1 / numpy.sqrt(diag))
            self.orbsym = symm.label_orb_symm(self.mol, irrep_name, self.mol.symm_orb, mo_coeff, s=s)

        if not hasattr(self.fcisolver, "orbsym") or not self.fcisolver.orbsym:
            ncore = self.ncore
            nocc = self.ncore + self.ncas
            self.fcisolver.orbsym = self.orbsym[ncore:nocc]
        logger.debug(self, "Active space irreps %s", str(self.fcisolver.orbsym))

        self.converged, self.e_tot, e_cas, self.ci, self.mo_coeff = _kern(
            self,
            mo_coeff,
            tol=self.conv_tol,
            conv_tol_grad=self.conv_tol_grad,
            macro=macro,
            micro=micro,
            ci0=ci0,
            callback=callback,
            verbose=self.verbose,
        )
        logger.note(self, "CASSCF energy = %.15g", self.e_tot)
        return self.e_tot, e_cas, self.ci, self.mo_coeff

    def get_occ(self, mo_energy, mo_coeff=None):
        ''' We cannot assume default mo_energy value, because the orbital
        energies are sorted after doing SCF.  But in this function, we need
        the orbital energies are grouped by symmetry irreps
        '''
        mol = self.mol
        nirrep = len(mol.symm_orb)
        if mo_coeff is not None:
            orbsym = symm.label_orb_symm(self, mol.irrep_id, mol.symm_orb,
                                         mo_coeff, self.get_ovlp(), False)
            orbsym = numpy.asarray(orbsym)
        else:
            orbsym = [numpy.repeat(ir, mol.symm_orb[i].shape[1])
                      for i, ir in enumerate(mol.irrep_id)]
            orbsym = numpy.hstack(orbsym)

        mo_occ = numpy.zeros_like(mo_energy)
        mo_e_left = []
        idx_e_left = []
        nelec_fix = 0
        for i, ir in enumerate(mol.irrep_id):
            irname = mol.irrep_name[i]
            ir_idx = numpy.where(orbsym == ir)[0]
            if irname in self.irrep_nelec:
                n = self.irrep_nelec[irname]
                e_idx = numpy.argsort(mo_energy[ir_idx])
                mo_occ[ir_idx[e_idx[:n//2]]] = 2
                nelec_fix += n
            else:
                idx_e_left.append(ir_idx)
        nelec_float = mol.nelectron - nelec_fix
        assert(nelec_float >= 0)
        if nelec_float > 0:
            idx_e_left = numpy.hstack(idx_e_left)
            mo_e_left = mo_energy[idx_e_left]
            mo_e_sort = numpy.argsort(mo_e_left)
            occ_idx = idx_e_left[mo_e_sort[:(nelec_float//2)]]
            mo_occ[occ_idx] = 2

        viridx = (mo_occ==0)
        if self.verbose < logger.INFO or viridx.sum() == 0:
            return mo_occ
        ehomo = max(mo_energy[mo_occ>0 ])
        elumo = min(mo_energy[mo_occ==0])
        noccs = []
        for i, ir in enumerate(mol.irrep_id):
            irname = mol.irrep_name[i]
            ir_idx = (orbsym == ir)

            noccs.append(int(mo_occ[ir_idx].sum()))
            if ehomo in mo_energy[ir_idx]:
                irhomo = irname
            if elumo in mo_energy[ir_idx]:
                irlumo = irname
        logger.info(self, 'HOMO (%s) = %.15g  LUMO (%s) = %.15g',
                    irhomo, ehomo, irlumo, elumo)
        if self.verbose >= logger.DEBUG:
            logger.debug(self, 'irrep_nelec = %s', noccs)
            _dump_mo_energy(mol, mo_energy, mo_occ, ehomo, elumo, orbsym)
        return mo_occ

def orbital_coeff(mol, fout, mo_coeff, spin='Alpha', symm=None, ene=None,
                  occ=None, ignore_h=False):
    from pyscf.symm import label_orb_symm
    if ignore_h:
        mol, mo_coeff = remove_high_l(mol, mo_coeff)
    aoidx = order_ao_index(mol)
    nmo = mo_coeff.shape[1]
    if symm is None:
        symm = ['A']*nmo
        if mol.symmetry:
            try:
                symm = label_orb_symm(mol, mol.irrep_name, mol.symm_orb,
                                      mo_coeff, tol=1e-5)
            except ValueError as e:
                logger.warn(mol, str(e))
    if ene is None:
        ene = numpy.arange(nmo)
    assert(spin == 'Alpha' or spin == 'Beta')
    if occ is None:
        occ = numpy.zeros(nmo)
        neleca, nelecb = mol.nelec
        if spin == 'Alpha':
            occ[:neleca] = 1
        else:
            occ[:nelecb] = 1
    fout.write('[MO]\n')
    for imo in range(nmo):
        fout.write(' Sym= %s\n' % symm[imo])
        fout.write(' Ene= %15.10g\n' % ene[imo])
        fout.write(' Spin= %s\n' % spin)
        fout.write(' Occup= %10.5f\n' % occ[imo])
        for i,j in enumerate(aoidx):
            fout.write(' %3d    %18.14g\n' % (i+1, mo_coeff[j,imo]))

def get_irrep_nelec(mol, mo_coeff, mo_occ, s=None):
    '''Electron numbers for each irreducible representation.

    Args:
        mol : an instance of :class:`Mole`
            To provide irrep_id, and spin-adapted basis
        mo_coeff : 2D ndarray
            Regular orbital coefficients, without grouping for irreps
        mo_occ : 1D ndarray
            Regular occupancy, without grouping for irreps

    Returns:
        irrep_nelec : dict
            The number of electrons for each irrep {'ir_name':int,...}.

    Examples:

    >>> mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='ccpvdz', symmetry=True, verbose=0)
    >>> mf = scf.RHF(mol)
    >>> mf.scf()
    -76.016789472074251
    >>> scf.hf_symm.get_irrep_nelec(mol, mf.mo_coeff, mf.mo_occ)
    {'A1': 6, 'A2': 0, 'B1': 2, 'B2': 2}
    '''
    orbsym = symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb, mo_coeff, s)
    orbsym = numpy.array(orbsym)
    irrep_nelec = dict([(mol.irrep_name[k], int(sum(mo_occ[orbsym==ir])))
                        for k, ir in enumerate(mol.irrep_id)])
    return irrep_nelec

    def test_kernel_symm(self):
        mol = gto.Mole()
        mol.verbose = 0
        mol.output = None
        mol.atom = [
            ['O', ( 0., 0.    , 0.   )],
            ['H', ( 0., -0.757, 0.587)],
            ['H', ( 0., 0.757 , 0.587)],]
        mol.basis = 'sto-3g'
        mol.symmetry = 1
        mol.build()
        m = scf.RHF(mol).run()

        norb = m.mo_coeff.shape[1]
        nelec = mol.nelectron - 2
        h1e = reduce(numpy.dot, (m.mo_coeff.T, scf.hf.get_hcore(mol), m.mo_coeff))
        eri = ao2mo.incore.full(m._eri, m.mo_coeff)
        orbsym = symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb, m.mo_coeff)

        myci = selected_ci_symm.SCI().set(orbsym=orbsym, select_cutoff=.5e-3)
        e1, c1 = myci.kernel(h1e, eri, norb, nelec)
        myci = direct_spin1_symm.FCISolver().set(orbsym=orbsym)
        e2, c2 = myci.kernel(h1e, eri, norb, nelec)

        self.assertAlmostEqual(e1, e2, 6)
        c2_cut = selected_ci.from_fci(c2, c1._strs, norb, nelec)
        self.assertAlmostEqual(abs(numpy.dot(c1.ravel(), c2_cut.ravel())), 1, 6)

        myci = selected_ci_symm.SCI().set(orbsym=orbsym, select_cutoff=1e-5)
        e1, c1 = myci.kernel(h1e, eri, norb, nelec, ci0=c2_cut, wfnsym=0)
        self.assertAlmostEqual(e1, e2, 9)
        c2_cut = selected_ci.from_fci(c2, c1._strs, norb, nelec)
        self.assertAlmostEqual(abs(numpy.dot(c1.ravel(), c2_cut.ravel())), 1, 9)

    def get_grad(self, mo_coeff, mo_occ, fock=None):
        mol = self.mol
        if not self.mol.symmetry:
            return rohf.ROHF.get_grad(self, mo_coeff, mo_occ, fock)

        if fock is None:
            dm1 = self.make_rdm1(mo_coeff, mo_occ)
            fock = self.get_hcore(mol) + self.get_veff(self.mol, dm1)
        orbsym = symm.label_orb_symm(self, mol.irrep_id, mol.symm_orb,
                                     mo_coeff, self.get_ovlp(), False)
        orbsym = numpy.asarray(orbsym)
        sym_allow = orbsym.reshape(-1,1) == orbsym

        occidxa = mo_occ > 0
        occidxb = mo_occ == 2
        viridxa = ~occidxa
        viridxb = ~occidxb
        focka = reduce(numpy.dot, (mo_coeff.T, fock[0], mo_coeff))
        fockb = reduce(numpy.dot, (mo_coeff.T, fock[1], mo_coeff))
        uniq_var_a = viridxa.reshape(-1,1) & occidxa
        uniq_var_b = viridxb.reshape(-1,1) & occidxb
        g = numpy.zeros_like(focka)
        g[uniq_var_a]  = focka[uniq_var_a]
        g[uniq_var_b] += fockb[uniq_var_b]
        g[~sym_allow] = 0
        return g[uniq_var_a | uniq_var_b].ravel()

def from_chkfile(filename, chkfile, tol=TOL, float_format=DEFAULT_FLOAT_FORMAT):
    '''Read SCF results from PySCF chkfile and transform 1-electron,
    2-electron integrals using the SCF orbitals.  The transformed integrals is
    written to FCIDUMP'''
    from pyscf import scf, symm
    with open(filename, 'w') as fout:
        mol, scf_rec = scf.chkfile.load_scf(chkfile)
        mo_coeff = numpy.array(scf_rec['mo_coeff'])
        nmo = mo_coeff.shape[1]
        if mol.symmetry:
            orbsym = symm.label_orb_symm(mol, mol.irrep_id,
                                         mol.symm_orb, mo_coeff, check=False)
            write_head(fout, nmo, mol.nelectron, mol.spin, orbsym)
        else:
            write_head(fout, nmo, mol.nelectron, mol.spin)

        eri = ao2mo.full(mol, mo_coeff, verbose=0)
        write_eri(fout, ao2mo.restore(8, eri, nmo), nmo, tol, float_format)

        t = mol.intor_symmetric('int1e_kin')
        v = mol.intor_symmetric('int1e_nuc')
        h = reduce(numpy.dot, (mo_coeff.T, t+v, mo_coeff))
        write_hcore(fout, h, nmo, tol, float_format)
        output_format = ' ' + float_format + '  0  0  0  0\n'
        fout.write(output_format % mol.energy_nuc())

    def analyze(self, verbose=logger.DEBUG):
        from pyscf.tools import dump_mat
        mo_energy = self.mo_energy
        mo_occ = self.mo_occ
        mo_coeff = self.mo_coeff
        log = logger.Logger(self.stdout, verbose)
        mol = self.mol
        nirrep = len(mol.irrep_id)
        ovlp_ao = self.get_ovlp()
        orbsym = symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb, mo_coeff,
                                     s=ovlp_ao)
        orbsym = numpy.array(orbsym)
        wfnsym = 0
        ndoccs = []
        nsoccs = []
        for k,ir in enumerate(mol.irrep_id):
            ndoccs.append(sum(orbsym[mo_occ==2] == ir))
            nsoccs.append(sum(orbsym[mo_occ==1] == ir))
            if nsoccs[k] % 2:
                wfnsym ^= ir
        if mol.groupname in ('Dooh', 'Coov'):
            log.info('TODO: total symmetry for %s', mol.groupname)
        else:
            log.info('total symmetry = %s',
                     symm.irrep_id2name(mol.groupname, wfnsym))
        log.info('occupancy for each irrep:  ' + (' %4s'*nirrep),
                 *mol.irrep_name)
        log.info('double occ                 ' + (' %4d'*nirrep), *ndoccs)
        log.info('single occ                 ' + (' %4d'*nirrep), *nsoccs)
        log.info('**** MO energy ****')
        irname_full = {}
        for k,ir in enumerate(mol.irrep_id):
            irname_full[ir] = mol.irrep_name[k]
        irorbcnt = {}
        for k, j in enumerate(orbsym):
            if j in irorbcnt:
                irorbcnt[j] += 1
            else:
                irorbcnt[j] = 1
            log.info('MO #%d (%s #%d), energy= %.15g occ= %g',
                     k+1, irname_full[j], irorbcnt[j], mo_energy[k], mo_occ[k])

        if verbose >= logger.DEBUG:
            label = mol.spheric_labels(True)
            molabel = []
            irorbcnt = {}
            for k, j in enumerate(orbsym):
                if j in irorbcnt:
                    irorbcnt[j] += 1
                else:
                    irorbcnt[j] = 1
                molabel.append('#%-d(%s #%d)' % (k+1, irname_full[j], irorbcnt[j]))
            log.debug(' ** MO coefficients **')
            dump_mat.dump_rec(mol.stdout, mo_coeff, label, molabel, start=1)

        dm = self.make_rdm1(mo_coeff, mo_occ)
        return self.mulliken_meta(mol, dm, s=ovlp_ao, verbose=verbose)

 def test_symmetrize_space(self):
     from pyscf import gto, symm, scf
     mol = gto.M(atom = 'C  0  0  0; H  1  1  1; H -1 -1  1; H  1 -1 -1; H -1  1 -1',
                 basis = 'sto3g', verbose=0)
     mf = scf.RHF(mol).run()
     mol.build(0, 0, symmetry='D2')
     mo = symm.symmetrize_space(mol, mf.mo_coeff)
     irreps = symm.label_orb_symm(mol, mol.irrep_name, mol.symm_orb, mo)
     self.assertEqual(irreps, ['A','A','A','B1','B1','B2','B2','B3','B3'])

 def __init__(self):
     self.__dict__.update(cis.__dict__)
     self.h1e = reduce(numpy.dot, (mo.T, scf.hf.get_hcore(mol), mo))
     self.eri = ao2mo.outcore.full_iofree(mol, mo)
     self.eci = None
     self.ci = None
     if mol.symmetry:
         self.orbsym = symm.label_orb_symm(mol, mol.irrep_id,
                                           mol.symm_orb, mo)
     self._keys = set(self.__dict__.keys())

def analyze(mf, verbose=logger.DEBUG, **kwargs):
    '''Analyze the given SCF object:  print orbital energies, occupancies;
    print orbital coefficients; Occupancy for each irreps; Mulliken population analysis
    '''
    from pyscf.lo import orth
    from pyscf.tools import dump_mat
    mol = mf.mol
    if not mol.symmetry:
        return hf.analyze(mf, verbose, **kwargs)

    mo_energy = mf.mo_energy
    mo_occ = mf.mo_occ
    mo_coeff = mf.mo_coeff
    log = logger.Logger(mf.stdout, verbose)
    nirrep = len(mol.irrep_id)
    ovlp_ao = mf.get_ovlp()
    orbsym = symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb, mo_coeff,
                                 s=ovlp_ao, check=False)
    orbsym = numpy.array(orbsym)
    wfnsym = 0
    noccs = [sum(orbsym[mo_occ>0]==ir) for ir in mol.irrep_id]
    log.note('total symmetry = %s', symm.irrep_id2name(mol.groupname, wfnsym))
    log.note('occupancy for each irrep:  ' + (' %4s'*nirrep), *mol.irrep_name)
    log.note('double occ                 ' + (' %4d'*nirrep), *noccs)
    log.note('**** MO energy ****')
    irname_full = {}
    for k,ir in enumerate(mol.irrep_id):
        irname_full[ir] = mol.irrep_name[k]
    irorbcnt = {}
    for k, j in enumerate(orbsym):
        if j in irorbcnt:
            irorbcnt[j] += 1
        else:
            irorbcnt[j] = 1
        log.note('MO #%d (%s #%d), energy= %.15g occ= %g',
                 k+1, irname_full[j], irorbcnt[j], mo_energy[k], mo_occ[k])

    if verbose >= logger.DEBUG:
        label = mol.spheric_labels(True)
        molabel = []
        irorbcnt = {}
        for k, j in enumerate(orbsym):
            if j in irorbcnt:
                irorbcnt[j] += 1
            else:
                irorbcnt[j] = 1
            molabel.append('#%-d(%s #%d)' % (k+1, irname_full[j], irorbcnt[j]))
        log.debug(' ** MO coefficients (expansion on meta-Lowdin AOs) **')
        orth_coeff = orth.orth_ao(mol, 'meta_lowdin', s=ovlp_ao)
        c = reduce(numpy.dot, (orth_coeff.T, ovlp_ao, mo_coeff))
        dump_mat.dump_rec(mf.stdout, c, label, molabel, start=1, **kwargs)

    dm = mf.make_rdm1(mo_coeff, mo_occ)
    return mf.mulliken_meta(mol, dm, s=ovlp_ao, verbose=log)

def label_symmetry_(mc, mo_coeff):
    #irrep_name = mc.mol.irrep_name
    irrep_name = mc.mol.irrep_id
    s = mc._scf.get_ovlp()
    try:
        mc.orbsym = symm.label_orb_symm(mc.mol, irrep_name,
                                        mc.mol.symm_orb, mo_coeff, s=s)
    except ValueError:
        logger.warn(mc, 'mc1step_symm symmetrizes input orbitals')
        ncore = mc.ncore
        nocc = mc.ncore + mc.ncas
        mo_cor = symm.symmetrize_space(mc.mol, mo_coeff[:,    :ncore], s=s)
        mo_act = symm.symmetrize_space(mc.mol, mo_coeff[:,ncore:nocc], s=s)
        mo_vir = symm.symmetrize_space(mc.mol, mo_coeff[:,nocc:     ], s=s)
        mo_coeff = numpy.hstack((mo_cor,mo_act,mo_vir))
        mc.orbsym = symm.label_orb_symm(mc.mol, irrep_name,
                                        mc.mol.symm_orb, mo_coeff, s=s)

    if not hasattr(mc.fcisolver, 'orbsym') or not mc.fcisolver.orbsym:
        ncore = mc.ncore
        nocc = mc.ncore + mc.ncas
        mc.fcisolver.orbsym = mc.orbsym[ncore:nocc]
    logger.debug(mc, 'Active space irreps %s', str(mc.fcisolver.orbsym))

    def get_grad(self, mo_coeff, mo_occ, fock=None):
        mol = self.mol
        if fock is None:
            dm1 = self.make_rdm1(mo_coeff, mo_occ)
            fock = self.get_hcore(mol) + self.get_veff(self.mol, dm1)
        orbsym = symm.label_orb_symm(self, mol.irrep_id, mol.symm_orb,
                                     mo_coeff, self.get_ovlp(), False)
        orbsym = numpy.asarray(orbsym)

        occidx = mo_occ > 0
        viridx = ~occidx
        g = reduce(numpy.dot, (mo_coeff[:,viridx].T.conj(), fock,
                               mo_coeff[:,occidx])) * 2
        sym_allow = orbsym[viridx].reshape(-1,1) == orbsym[occidx]
        g[~sym_allow] = 0
        return g.ravel()

def write_integrals(xci, orb):
    mol = xci.mol
    orbsym = symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb, orb)
    h1e = reduce(numpy.dot, (orb.T, xci.get_hcore(), orb))
    norb = orb.shape[1]
    if xci._eri is not None:
        h2e = ao2mo.restore(1, ao2mo.full(xci._eri, orb), norb)
    else:
        h2e = ao2mo.restore(1, ao2mo.full(mol, orb), norb)

    with h5py.File(xci.integralfile, 'w') as f:
        f['h1e']    = h1e
        f['h2e']    = h2e
        f['norb'  ] = numpy.array(norb, dtype=numpy.int32)
        f['group' ] = mol.groupname
        f['orbsym'] = numpy.asarray(orbsym, dtype=numpy.int32)
        f['ecore' ] = mol.energy_nuc()

def get_orbsym(mol, mo_coeff, s=None, check=False):
    if mo_coeff is None:
        orbsym = numpy.hstack([[ir] * mol.symm_orb[i].shape[1]
                               for i, ir in enumerate(mol.irrep_id)])
    elif getattr(mo_coeff, 'orbsym', None) is not None:
        orbsym = mo_coeff.orbsym
    else:
        nao = mo_coeff.shape[0] // 2
        if isinstance(s, numpy.ndarray):
            assert(s.size == nao**2 or numpy.allclose(s[:nao,:nao], s[nao:,nao:]))
            s = s[:nao,:nao]
        mo_a = mo_coeff[:nao].copy()
        mo_b = mo_coeff[nao:]
        zero_alpha_idx = numpy.linalg.norm(mo_a, axis=0) < 1e-7
        mo_a[:,zero_alpha_idx] = mo_b[:,zero_alpha_idx]
        orbsym = symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb,
                                     mo_a, s, check)
    return numpy.asarray(orbsym)

def label_symmetry_(mc, mo_coeff, ci0=None):
    log = logger.Logger(mc.stdout, mc.verbose)
    #irrep_name = mc.mol.irrep_name
    irrep_name = mc.mol.irrep_id
    s = mc._scf.get_ovlp()
    ncore = mc.ncore
    nocc = ncore + mc.ncas
    try:
        orbsym = scf.hf_symm.get_orbsym(mc._scf.mol, mo_coeff, s, True)
    except ValueError:
        log.warn('mc1step_symm symmetrizes input orbitals')
        mo_cor = symm.symmetrize_space(mc.mol, mo_coeff[:,    :ncore], s=s, check=False)
        mo_act = symm.symmetrize_space(mc.mol, mo_coeff[:,ncore:nocc], s=s, check=False)
        mo_vir = symm.symmetrize_space(mc.mol, mo_coeff[:,nocc:     ], s=s, check=False)
        mo_coeff = numpy.hstack((mo_cor,mo_act,mo_vir))
        orbsym = symm.label_orb_symm(mc.mol, irrep_name,
                                     mc.mol.symm_orb, mo_coeff, s=s)
    mo_coeff_with_orbsym = lib.tag_array(mo_coeff, orbsym=orbsym)

    active_orbsym = getattr(mc.fcisolver, 'orbsym', [])
    if (not getattr(active_orbsym, '__len__', None)) or len(active_orbsym) == 0:
        mc.fcisolver.orbsym = orbsym[ncore:nocc]
    log.debug('Active space irreps %s', str(mc.fcisolver.orbsym))

    wfnsym = 0
    if getattr(mc.fcisolver, 'wfnsym', None) is not None:
        wfnsym = mc.fcisolver.wfnsym

    elif ci0 is None:
        # Guess wfnsym based on HF determinant.  mo_coeff may not be HF
        # canonical orbitals.  Some checks are needed to ensure that mo_coeff
        # are derived from the symmetry adapted SCF calculations.
        if mo_coeff is mc._scf.mo_coeff:
            wfnsym = 0
            for ir in orbsym[mc._scf.mo_occ == 1]:
                wfnsym ^= ir
            mc.fcisolver.wfnsym = wfnsym
            log.debug('Set CASCI wfnsym %s based on HF determinant', wfnsym)
        elif getattr(mo_coeff, 'orbsym', None) is not None:  # It may be reordered SCF orbitals
            cas_orb = mo_coeff[:,ncore:nocc]
            s = reduce(numpy.dot, (cas_orb.T, mc._scf.get_ovlp(), mc._scf.mo_coeff))
            if numpy.all(numpy.max(s, axis=1) > 1-1e-9):
                idx = numpy.argmax(s, axis=1)
                cas_orbsym = orbsym[ncore:nocc]
                cas_occ = mc._scf.mo_occ[idx]
                wfnsym = 0
                for ir in cas_orbsym[cas_occ == 1]:
                    wfnsym ^= ir
                mc.fcisolver.wfnsym = wfnsym
                log.debug('Active space are constructed from canonical SCF '
                          'orbitals %s', idx)
                log.debug('Set CASCI wfnsym %s based on HF determinant', wfnsym)

    elif getattr(mc.fcisolver, 'guess_wfnsym', None):
        wfnsym = mc.fcisolver.guess_wfnsym(mc.ncas, mc.nelecas, ci0, verbose=log)
        log.debug('CASCI wfnsym %s (based on CI initial guess)', wfnsym)

    if isinstance(wfnsym, (int, numpy.integer)):
        wfnsym = symm.irrep_id2name(mc.mol.groupname, wfnsym)

    log.info('Active space CI wfn symmetry = %s', wfnsym)

    return mo_coeff_with_orbsym

    def analyze(self, verbose=logger.DEBUG):
        from pyscf.lo import orth
        from pyscf.tools import dump_mat
        mo_energy = self.mo_energy
        mo_occ = self.mo_occ
        mo_coeff = self.mo_coeff
        log = logger.Logger(self.stdout, verbose)
        mol = self.mol
        nirrep = len(mol.irrep_id)
        ovlp_ao = self.get_ovlp()
        orbsym = symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb, mo_coeff,
                                     s=ovlp_ao)
        orbsym = numpy.array(orbsym)
        wfnsym = 0
        ndoccs = []
        nsoccs = []
        for k,ir in enumerate(mol.irrep_id):
            ndoccs.append(sum(orbsym[mo_occ==2] == ir))
            nsoccs.append(sum(orbsym[mo_occ==1] == ir))
            if nsoccs[k] % 2:
                wfnsym ^= ir
        if mol.groupname in ('Dooh', 'Coov'):
            log.info('TODO: total symmetry for %s', mol.groupname)
        else:
            log.info('total symmetry = %s',
                     symm.irrep_id2name(mol.groupname, wfnsym))
        log.info('occupancy for each irrep:  ' + (' %4s'*nirrep),
                 *mol.irrep_name)
        log.info('double occ                 ' + (' %4d'*nirrep), *ndoccs)
        log.info('single occ                 ' + (' %4d'*nirrep), *nsoccs)
        log.info('**** MO energy ****')
        irname_full = {}
        for k,ir in enumerate(mol.irrep_id):
            irname_full[ir] = mol.irrep_name[k]
        irorbcnt = {}
        if self._focka_ao is None:
            for k, j in enumerate(orbsym):
                if j in irorbcnt:
                    irorbcnt[j] += 1
                else:
                    irorbcnt[j] = 1
                log.note('MO #%-3d (%s #%-2d), energy= %-18.15g occ= %g',
                         k+1, irname_full[j], irorbcnt[j],
                         mo_energy[k], mo_occ[k])
        else:
            mo_ea = numpy.einsum('ik,ik->k', mo_coeff, self._focka_ao.dot(mo_coeff))
            mo_eb = numpy.einsum('ik,ik->k', mo_coeff, self._fockb_ao.dot(mo_coeff))
            log.note('                          Roothaan           | alpha              | beta')
            for k, j in enumerate(orbsym):
                if j in irorbcnt:
                    irorbcnt[j] += 1
                else:
                    irorbcnt[j] = 1
                log.note('MO #%-4d(%-3s #%-2d) energy= %-18.15g | %-18.15g | %-18.15g occ= %g',
                         k+1, irname_full[j], irorbcnt[j],
                         mo_energy[k], mo_ea[k], mo_eb[k], mo_occ[k])

        if verbose >= logger.DEBUG:
            label = mol.spheric_labels(True)
            molabel = []
            irorbcnt = {}
            for k, j in enumerate(orbsym):
                if j in irorbcnt:
                    irorbcnt[j] += 1
                else:
                    irorbcnt[j] = 1
                molabel.append('#%-d(%s #%d)' % (k+1, irname_full[j], irorbcnt[j]))
            log.debug(' ** MO coefficients (expansion on meta-Lowdin AOs) **')
            orth_coeff = orth.orth_ao(mol, 'meta_lowdin', s=ovlp_ao)
            c = reduce(numpy.dot, (orth_coeff.T, ovlp_ao, mo_coeff))
            dump_mat.dump_rec(self.stdout, c, label, molabel, start=1)

        dm = self.make_rdm1(mo_coeff, mo_occ)
        return self.mulliken_meta(mol, dm, s=ovlp_ao, verbose=verbose)