def analyze(mf, verbose=logger.DEBUG, **kwargs):
    '''Analyze the given SCF object:  print orbital energies, occupancies;
    print orbital coefficients; Mulliken population analysis
    '''
    from pyscf.lo import orth
    from pyscf.tools import dump_mat
    mo_energy = mf.mo_energy
    mo_occ = mf.mo_occ
    mo_coeff = mf.mo_coeff
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(mf.stdout, verbose)

    log.note('**** MO energy ****')
    if mf._focka_ao is None:
        for i,c in enumerate(mo_occ):
            log.note('MO #%-3d energy= %-18.15g occ= %g', i+1, mo_energy[i], c)
    else:
        mo_ea = numpy.einsum('ik,ik->k', mo_coeff, mf._focka_ao.dot(mo_coeff))
        mo_eb = numpy.einsum('ik,ik->k', mo_coeff, mf._fockb_ao.dot(mo_coeff))
        log.note('                Roothaan           | alpha              | beta')
        for i,c in enumerate(mo_occ):
            log.note('MO #%-3d energy= %-18.15g | %-18.15g | %-18.15g occ= %g',
                     i+1, mo_energy[i], mo_ea[i], mo_eb[i], c)
    ovlp_ao = mf.get_ovlp()
    if verbose >= logger.DEBUG:
        log.debug(' ** MO coefficients (expansion on meta-Lowdin AOs) **')
        label = mf.mol.spheric_labels(True)
        orth_coeff = orth.orth_ao(mf.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, start=1, **kwargs)
    dm = mf.make_rdm1(mo_coeff, mo_occ)
    return mf.mulliken_meta(mf.mol, dm, s=s, verbose=log)

def analyze(mf, verbose=logger.DEBUG):
    '''Analyze the given SCF object:  print orbital energies, occupancies;
    print orbital coefficients; Mulliken population analysis; Diople moment.
    '''
    from pyscf.lo import orth
    from pyscf.tools import dump_mat
    mo_energy = mf.mo_energy
    mo_occ = mf.mo_occ
    mo_coeff = mf.mo_coeff
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(mf.stdout, verbose)

    log.note('**** MO energy ****')
    for i,c in enumerate(mo_occ):
        log.note('MO #%-3d energy= %-18.15g occ= %g', i+1, mo_energy[i], c)
    ovlp_ao = mf.get_ovlp()
    if verbose >= logger.DEBUG:
        log.debug(' ** MO coefficients (expansion on meta-Lowdin AOs) **')
        label = mf.mol.spheric_labels(True)
        orth_coeff = orth.orth_ao(mf.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, start=1)
    dm = mf.make_rdm1(mo_coeff, mo_occ)
    return (mf.mulliken_meta(mf.mol, dm, s=ovlp_ao, verbose=log),
            mf.dip_moment(mf.mol, dm, verbose=log))

def analyze(mf, verbose=logger.DEBUG):
    """Analyze the given SCF object:  print orbital energies, occupancies;
    print orbital coefficients; Mulliken population analysis
    """
    from pyscf.tools import dump_mat

    mo_energy = mf.mo_energy
    mo_occ = mf.mo_occ
    mo_coeff = mf.mo_coeff
    log = logger.Logger(mf.stdout, verbose)
    ss, s = mf.spin_square((mo_coeff[0][:, mo_occ[0] > 0], mo_coeff[1][:, mo_occ[1] > 0]), mf.get_ovlp())
    log.info("multiplicity <S^2> = %.8g  2S+1 = %.8g", ss, s)

    log.info("**** MO energy ****")
    for i in range(mo_energy[0].__len__()):
        if mo_occ[0][i] > 0:
            log.info("alpha occupied MO #%d energy = %.15g occ= %g", i + 1, mo_energy[0][i], mo_occ[0][i])
        else:
            log.info("alpha virtual MO #%d energy = %.15g occ= %g", i + 1, mo_energy[0][i], mo_occ[0][i])
    for i in range(mo_energy[1].__len__()):
        if mo_occ[1][i] > 0:
            log.info("beta occupied MO #%d energy = %.15g occ= %g", i + 1, mo_energy[1][i], mo_occ[1][i])
        else:
            log.info("beta virtual MO #%d energy = %.15g occ= %g", i + 1, mo_energy[1][i], mo_occ[1][i])
    if mf.verbose >= logger.DEBUG:
        log.debug(" ** MO coefficients for alpha spin **")
        label = mf.mol.spheric_labels(True)
        dump_mat.dump_rec(mf.stdout, mo_coeff[0], label, start=1)
        log.debug(" ** MO coefficients for beta spin **")
        dump_mat.dump_rec(mf.stdout, mo_coeff[1], label, start=1)

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

def analyze(mf, verbose=logger.DEBUG, **kwargs):
    from pyscf.tools import dump_mat
    log = logger.new_logger(mf, verbose)
    mo_energy = mf.mo_energy
    mo_occ = mf.mo_occ
    mo_coeff = mf.mo_coeff

    log.info('**** MO energy ****')
    for i in range(len(mo_energy)):
        if mo_occ[i] > 0:
            log.info('occupied MO #%d energy= %.15g occ= %g', \
                     i+1, mo_energy[i], mo_occ[i])
        else:
            log.info('virtual MO #%d energy= %.15g occ= %g', \
                     i+1, mo_energy[i], mo_occ[i])
    mol = mf.mol
    if mf.verbose >= logger.DEBUG1:
        log.debug(' ** MO coefficients of large component of postive state (real part) **')
        label = mol.spinor_labels()
        n2c = mo_coeff.shape[0] // 2
        dump_mat.dump_rec(mf.stdout, mo_coeff[n2c:,:n2c].real, label, start=1)
    dm = mf.make_rdm1(mo_coeff, mo_occ)
    pop_chg = mf.mulliken_pop(mol, dm, mf.get_ovlp(), log)
    dip = mf.dip_moment(mol, dm, verbose=log)
    return pop_chg, dip

def analyze(mf, verbose=logger.DEBUG):
    '''Analyze the given SCF object:  print orbital energies, occupancies;
    print orbital coefficients; Mulliken population analysis
    '''
    from pyscf.tools import dump_mat
    mo_energy = mf.mo_energy
    mo_occ = mf.mo_occ
    mo_coeff = mf.mo_coeff
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(mf.stdout, verbose)

    log.info('**** MO energy ****')
    for i in range(len(mo_energy)):
        if mo_occ[i] > 0:
            log.info('occupied MO #%d energy= %.15g occ= %g',
                     i+1, mo_energy[i], mo_occ[i])
        else:
            log.info('virtual MO #%d energy= %.15g occ= %g',
                     i+1, mo_energy[i], mo_occ[i])
    if verbose >= logger.DEBUG:
        log.debug(' ** MO coefficients **')
        label = ['%d%3s %s%-4s' % x for x in mf.mol.spheric_labels()]
        dump_mat.dump_rec(mf.stdout, mo_coeff, label, start=1)
    dm = mf.make_rdm1(mo_coeff, mo_occ)
    return mf.mulliken_pop(mf.mol, dm, mf.get_ovlp(), log)

    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 analyze(casscf, mo_coeff=None, ci=None, verbose=logger.INFO):
    from pyscf.tools import dump_mat
    from pyscf.mcscf import addons
    if mo_coeff is None: mo_coeff = casscf.mo_coeff
    if ci is None: ci = casscf.ci
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(casscf.stdout, verbose)
    nelecas = casscf.nelecas
    ncas = casscf.ncas
    ncore = casscf.ncore
    nocc = ncore + ncas

    casdm1a, casdm1b = casscf.fcisolver.make_rdm1s(ci, ncas, nelecas)
    mocore = mo_coeff[:,:ncore]
    mocas = mo_coeff[:,ncore:nocc]
    dm1b = numpy.dot(mocore, mocore.T)
    dm1a = dm1b + reduce(numpy.dot, (mocas, casdm1a, mocas.T))
    dm1b += reduce(numpy.dot, (mocas, casdm1b, mocas.T))

    if log.verbose >= logger.INFO:
        label = ['%d%3s %s%-4s' % x for x in casscf.mol.spheric_labels()]
        if log.verbose >= logger.DEBUG:
            log.info('alpha density matrix (on AO)')
            dump_mat.dump_tri(log.stdout, dm1a, label)
            log.info('beta density matrix (on AO)')
            dump_mat.dump_tri(log.stdout, dm1b, label)

        # note the last two args of ._eig for mc1step_symm
        occ, ucas = casscf._eig(-(casdm1a+casdm1b), ncore, nocc)
        log.info('Natural occ %s', str(-occ))
        for i, k in enumerate(numpy.argmax(abs(ucas), axis=0)):
            if ucas[k,i] < 0:
                ucas[:,i] *= -1
        mo_cas = numpy.dot(mo_coeff[:,ncore:nocc], ucas)
        log.info('Natural orbital in CAS space')
        dump_mat.dump_rec(log.stdout, mo_cas, label, start=1)

        s = reduce(numpy.dot, (casscf.mo_coeff.T, casscf._scf.get_ovlp(),
                               casscf._scf.mo_coeff))
        idx = numpy.argwhere(abs(s)>.4)
        for i,j in idx:
            log.info('<mo-mcscf|mo-hf> %d, %d, %12.8f', i+1, j+1, s[i,j])

        log.info('** Largest CI components **')
        log.info(' string alpha, string beta, CI coefficients')
        for c,ia,ib in fci.addons.large_ci(ci, casscf.ncas, casscf.nelecas):
            log.info('  %9s    %9s    %.12f', ia, ib, c)

        dm1 = dm1a + dm1b
        s = casscf._scf.get_ovlp()
        casscf._scf.mulliken_pop(casscf.mol, dm1, s, verbose=log)
        casscf._scf.mulliken_pop_meta_lowdin_ao(casscf.mol, dm1, verbose=log)
    return dm1a, dm1b

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 analyze(mf, verbose=logger.DEBUG):
    '''Analyze the given SCF object:  print orbital energies, occupancies;
    print orbital coefficients; Occupancy for each irreps; Mulliken population analysis
    '''
    from pyscf.tools import dump_mat
    mo_energy = mf.mo_energy
    mo_occ = mf.mo_occ
    mo_coeff = mf.mo_coeff
    log = pyscf.lib.logger.Logger(mf.stdout, verbose)
    mol = mf.mol
    nirrep = len(mol.irrep_id)
    ovlp_ao = mf.get_ovlp()
    orbsym = pyscf.symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb,
                                       mo_coeff, s=ovlp_ao)
    orbsym = numpy.array(orbsym)
    tot_sym = 0
    noccs = [sum(orbsym[mo_occ>0]==ir) for ir in mol.irrep_id]
    log.info('total symmetry = %s',
             pyscf.symm.irrep_name(mol.groupname, tot_sym))
    log.info('occupancy for each irrep:  ' + (' %4s'*nirrep), *mol.irrep_name)
    log.info('double occ                 ' + (' %4d'*nirrep), *noccs)
    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 = ['%d%3s %s%-4s' % x for x in mol.spheric_labels()]
        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 = mf.make_rdm1(mo_coeff, mo_occ)
    return mf.mulliken_pop(mol, dm, ovlp_ao, log)

def analyze(mf, verbose=logger.DEBUG, with_meta_lowdin=WITH_META_LOWDIN,
            **kwargs):
    '''Analyze the given SCF object:  print orbital energies, occupancies;
    print orbital coefficients; Mulliken population analysis
    '''
    from pyscf.lo import orth
    from pyscf.tools import dump_mat
    mo_energy = mf.mo_energy
    mo_occ = mf.mo_occ
    mo_coeff = mf.mo_coeff
    log = logger.new_logger(mf, verbose)
    if log.verbose >= logger.NOTE:
        log.note('**** MO energy ****')
        if getattr(mo_energy, 'mo_ea', None) is not None:
            mo_ea = mo_energy.mo_ea
            mo_eb = mo_energy.mo_eb
            log.note('                Roothaan           | alpha              | beta')
            for i,c in enumerate(mo_occ):
                log.note('MO #%-3d energy= %-18.15g | %-18.15g | %-18.15g occ= %g',
                         i+MO_BASE, mo_energy[i], mo_ea[i], mo_eb[i], c)
        else:
            for i,c in enumerate(mo_occ):
                log.note('MO #%-3d energy= %-18.15g occ= %g',
                         i+MO_BASE, mo_energy[i], c)

    ovlp_ao = mf.get_ovlp()
    if log.verbose >= logger.DEBUG:
        label = mf.mol.ao_labels()
        if with_meta_lowdin:
            log.debug(' ** MO coefficients (expansion on meta-Lowdin AOs) **')
            orth_coeff = orth.orth_ao(mf.mol, 'meta_lowdin', s=ovlp_ao)
            c = reduce(numpy.dot, (orth_coeff.T, ovlp_ao, mo_coeff))
        else:
            log.debug(' ** MO coefficients (expansion on AOs) **')
            c = mo_coeff
        dump_mat.dump_rec(mf.stdout, c, label, start=MO_BASE, **kwargs)
    dm = mf.make_rdm1(mo_coeff, mo_occ)
    if with_meta_lowdin:
        pop_and_charge = mf.mulliken_meta(mf.mol, dm, s=ovlp_ao, verbose=log)
    else:
        pop_and_charge = mf.mulliken_pop(mf.mol, dm, s=ovlp_ao, verbose=log)
    dip = mf.dip_moment(mf.mol, dm, verbose=log)
    return pop_and_charge, dip

def analyze(mf, verbose=logger.DEBUG, **kwargs):
    '''Analyze the given SCF object:  print orbital energies, occupancies;
    print orbital coefficients; Mulliken population analysis; Dipole moment
    '''
    from pyscf.lo import orth
    from pyscf.tools import dump_mat
    mo_energy = mf.mo_energy
    mo_occ = mf.mo_occ
    mo_coeff = mf.mo_coeff
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(mf.stdout, verbose)
    ss, s = mf.spin_square((mo_coeff[0][:,mo_occ[0]>0],
                            mo_coeff[1][:,mo_occ[1]>0]), mf.get_ovlp())
    log.note('multiplicity <S^2> = %.8g  2S+1 = %.8g', ss, s)

    log.note('**** MO energy ****')
    log.note('                             alpha | beta                alpha | beta')
    for i in range(mo_occ.shape[1]):
        log.note('MO #%-3d energy= %-18.15g | %-18.15g occ= %g | %g',
                 i+1, mo_energy[0][i], mo_energy[1][i],
                 mo_occ[0][i], mo_occ[1][i])
    ovlp_ao = mf.get_ovlp()
    if verbose >= logger.DEBUG:
        log.debug(' ** MO coefficients (expansion on meta-Lowdin AOs) for alpha spin **')
        label = mf.mol.spheric_labels(True)
        orth_coeff = orth.orth_ao(mf.mol, 'meta_lowdin', s=ovlp_ao)
        c_inv = numpy.dot(orth_coeff.T, ovlp_ao)
        dump_mat.dump_rec(mf.stdout, c_inv.dot(mo_coeff[0]), label, start=1,
                          **kwargs)
        log.debug(' ** MO coefficients (expansion on meta-Lowdin AOs) for beta spin **')
        dump_mat.dump_rec(mf.stdout, c_inv.dot(mo_coeff[1]), label, start=1,
                          **kwargs)

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

def analyze(mf, verbose=logger.DEBUG, with_meta_lowdin=WITH_META_LOWDIN,
            **kwargs):
    '''Analyze the given SCF object:  print orbital energies, occupancies;
    print orbital coefficients; Mulliken population analysis; Dipole moment
    '''
    from pyscf.lo import orth
    from pyscf.tools import dump_mat
    mo_energy = mf.mo_energy
    mo_occ = mf.mo_occ
    mo_coeff = mf.mo_coeff
    nmo = len(mo_occ[0])
    log = logger.new_logger(mf, verbose)
    if log.verbose >= logger.NOTE:
        log.note('**** MO energy ****')
        log.note('                             alpha | beta                alpha | beta')
        for i in range(nmo):
            log.note('MO #%-3d energy= %-18.15g | %-18.15g occ= %g | %g',
                     i+MO_BASE, mo_energy[0][i], mo_energy[1][i],
                     mo_occ[0][i], mo_occ[1][i])

    ovlp_ao = mf.get_ovlp()
    if log.verbose >= logger.DEBUG:
        label = mf.mol.ao_labels()
        if with_meta_lowdin:
            log.debug(' ** MO coefficients (expansion on meta-Lowdin AOs) for alpha spin **')
            orth_coeff = orth.orth_ao(mf.mol, 'meta_lowdin', s=ovlp_ao)
            c_inv = numpy.dot(orth_coeff.T, ovlp_ao)
            dump_mat.dump_rec(mf.stdout, c_inv.dot(mo_coeff[0]), label,
                              start=MO_BASE, **kwargs)
            log.debug(' ** MO coefficients (expansion on meta-Lowdin AOs) for beta spin **')
            dump_mat.dump_rec(mf.stdout, c_inv.dot(mo_coeff[1]), label,
                              start=MO_BASE, **kwargs)
        else:
            log.debug(' ** MO coefficients (expansion on AOs) for alpha spin **')
            dump_mat.dump_rec(mf.stdout, mo_coeff[0], label,
                              start=MO_BASE, **kwargs)
            log.debug(' ** MO coefficients (expansion on AOs) for beta spin **')
            dump_mat.dump_rec(mf.stdout, mo_coeff[1], label,
                              start=MO_BASE, **kwargs)

    dm = mf.make_rdm1(mo_coeff, mo_occ)
    if with_meta_lowdin:
        return (mf.mulliken_meta(mf.mol, dm, s=ovlp_ao, verbose=log),
                mf.dip_moment(mf.mol, dm, verbose=log))
    else:
        return (mf.mulliken_pop(mf.mol, dm, s=ovlp_ao, verbose=log),
                mf.dip_moment(mf.mol, dm, verbose=log))

def cas_natorb(mc, mo_coeff=None, ci=None, eris=None, sort=False,
               casdm1=None, verbose=None):
    '''Transform active orbitals to natrual orbitals, and update the CI wfn

    Args:
        mc : a CASSCF/CASCI object or RHF object

    Kwargs:
        sort : bool
            Sort natural orbitals wrt the occupancy.  Be careful with this
            option since the resultant natural orbitals might have the
            different symmetry to the irreps indicated by CASSCF.orbsym

    Returns:
        A tuple, the first item is natural orbitals, the second is updated CI
        coefficients, the third is the natural occupancy associated to the
        natural orbitals.
    '''
    from pyscf.lo import orth
    from pyscf.tools import dump_mat
    from pyscf.tools.mo_mapping import mo_1to1map
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(mc.stdout, mc.verbose)
    if mo_coeff is None: mo_coeff = mc.mo_coeff
    if ci is None: ci = mc.ci
    ncore = mc.ncore
    ncas = mc.ncas
    nocc = ncore + ncas
    nelecas = mc.nelecas
    if casdm1 is None:
        casdm1 = mc.fcisolver.make_rdm1(ci, ncas, nelecas)
    # orbital symmetry is reserved in this _eig call
    occ, ucas = mc._eig(-casdm1, ncore, nocc)
    if sort:
        idx = numpy.argsort(occ)
        occ = occ[idx]
        ucas = ucas[:,idx]
# restore phase
# where_natorb gives the location of the natural orbital for the input cas
# orbitals.  gen_strings4orblist map thes sorted strings (on CAS orbital) to
# the unsorted determinant strings (on natural orbital). e.g.  (3o,2e) system
#       CAS orbital      1  2  3
#       natural orbital  3  1  2        <= by mo_1to1map
#       CASorb-strings   0b011, 0b101, 0b110
#                    ==  (1,2), (1,3), (2,3)
#       natorb-strings   (3,1), (3,2), (1,2)
#                    ==  0B101, 0B110, 0B011    <= by gen_strings4orblist
# then argsort to translate the string representation to the address
#       [2(=0B011), 0(=0B101), 1(=0B110)]
# to indicate which CASorb-strings address to be loaded in each natorb-strings slot
    where_natorb = mo_1to1map(ucas)
    for i, k in enumerate(where_natorb):
        if ucas[i,k] < 0:
            ucas[:,k] *= -1

    occ = -occ
    mo_occ = numpy.zeros(mo_coeff.shape[1])
    mo_occ[:ncore] = 2
    mo_occ[ncore:nocc] = occ

    if isinstance(ci, numpy.ndarray):
        fcivec = fci.addons.transform_ci_for_orbital_rotation(ci, ncas, nelecas, ucas)
    elif isinstance(ci, (tuple, list)) and isinstance(ci[0], numpy.ndarray):
        # for state-average eigenfunctions
        fcivec = [fci.addons.transform_ci_for_orbital_rotation(x, ncas, nelecas, ucas)
                  for x in ci]
    else:
        log.info('FCI vector not available, call CASCI for wavefunction')
        mocas = mo_coeff1[:,ncore:nocc]
        h1eff = reduce(numpy.dot, (mocas.T, mc.get_hcore(), mocas))
        if eris is not None and hasattr(eris, 'ppaa'):
            h1eff += reduce(numpy.dot, (ucas.T, eris.vhf_c[ncore:nocc,ncore:nocc], ucas))
            aaaa = ao2mo.restore(4, eris.ppaa[ncore:nocc,ncore:nocc,:,:], ncas)
            aaaa = ao2mo.incore.full(aaaa, ucas)
        else:
            dm_core = numpy.dot(mo_coeff[:,:ncore]*2, mo_coeff[:,:ncore].T)
            vj, vk = mc._scf.get_jk(mc.mol, dm_core)
            h1eff += reduce(numpy.dot, (mocas.T, vj-vk*.5, mocas))
            aaaa = ao2mo.kernel(mc.mol, mocas)
        max_memory = max(400, mc.max_memory-lib.current_memory()[0])
        e_cas, fcivec = mc.fcisolver.kernel(h1eff, aaaa, ncas, nelecas,
                                            max_memory=max_memory, verbose=log)
        log.debug('In Natural orbital, CI energy = %.12g', e_cas)

    mo_coeff1 = mo_coeff.copy()
    mo_coeff1[:,ncore:nocc] = numpy.dot(mo_coeff[:,ncore:nocc], ucas)
    if log.verbose >= logger.INFO:
        ovlp_ao = mc._scf.get_ovlp()
        log.debug('where_natorb %s', str(where_natorb))
        log.info('Natural occ %s', str(occ))
        log.info('Natural orbital (expansion on meta-Lowdin AOs) in CAS space')
        label = mc.mol.spheric_labels(True)
        orth_coeff = orth.orth_ao(mc.mol, 'meta_lowdin', s=ovlp_ao)
        mo_cas = reduce(numpy.dot, (orth_coeff.T, ovlp_ao, mo_coeff1[:,ncore:nocc]))
        dump_mat.dump_rec(log.stdout, mo_cas, label, start=1)

        if mc._scf.mo_coeff is not None:
            s = reduce(numpy.dot, (mo_coeff1[:,ncore:nocc].T,
#.........这里部分代码省略.........

#.........这里部分代码省略.........
    nao = dm.shape[0]
    nimp = len(baslst)
    log.debug('*** decompose density matrix')
    log.debug('orth AO method = %s', orth_method)
    log.debug('embedding AO list = %s', str(baslst))
    if orth_method is not None:
        corth = lo.orth.orth_ao(mol, method=orth_method, s=s)
        cinv = numpy.dot(corth.T, s)
        dm = reduce(numpy.dot, (cinv, dm, cinv.T))
    else:
        corth = numpy.eye(nao)

    baslst = numpy.asarray(baslst)
    notimp = numpy.asarray([i for i in range(nao) if i not in baslst])
    occi, ui = scipy.linalg.eigh(-dm[baslst[:,None],baslst])
    occi *= -1
    idxi = numpy.argsort(abs(occi-1))
    log.debug('entanglement weight occ = %s', str(occi[idxi]))
    occb, ub = scipy.linalg.eigh(dm[notimp[:,None],notimp])
    idxb = numpy.argsort(abs(occb-1))  # sort by entanglement
    occb = occb[idxb]
    ub = ub[:,idxb]

# guess ncas and nelecas
    nb = ((occb > occ_cutoff) & (occb < 2-occ_cutoff)).sum()
    log.debug('bath weight occ = %s', occb[:nb])
    cum_nelec = numpy.cumsum(occb[:nb]) + occi.sum()
    cum_nelec = numpy.append(occi.sum(), cum_nelec)
    log.debug('Active space cum nelec imp|[baths] = %f |%s',
              cum_nelec[0], cum_nelec[1:])
    ne_error = abs(cum_nelec.round() - cum_nelec)
    nb4cas = nb
    for i in range(nb):
        if (ne_error[i] < nelec_tol and
# whether all baths next to ith bath are less important
            (occb[i] < nelec_tol or occb[i] > 2-nelec_tol)):
            nb4cas = i
            break
    ncas = nb4cas + nimp
    nelecas = int(cum_nelec[nb4cas].round())
    ncore = (mol.nelectron - nelecas) // 2
    log.info('From DMET guess, ncas = %d  nelecas = %d  ncore = %d',
             ncas, nelecas, ncore)

    log.debug('DMET impurity and bath orbitals on orthogonal AOs')
    log.debug('DMET %d impurity sites/occ', nimp)
    if log.verbose >= logger.DEBUG1:
        label = mol.spherical_labels(True)
        occ_label = ['#%d/%.5f'%(i+1,x) for i,x in enumerate(occi)]
        #dump_mat.dump_rec(mol.stdout, numpy.dot(corth[:,baslst], ui),
        #                  label=label, label2=occ_label, start=1)
        dump_mat.dump_rec(mol.stdout, ui, label=[label[i] for i in baslst],
                          label2=occ_label, start=1)

    log.debug('DMET %d entangled baths/occ', nb)
    if log.verbose >= logger.DEBUG1:
        occ_label = ['#%d/%.5f'%(i+1,occb[i]) for i in range(nb)]
        #dump_mat.dump_rec(mol.stdout, numpy.dot(corth[:,notimp], ub[:,:nb]),
        #                  label=label, label2=occ_label, start=1)
        dump_mat.dump_rec(mol.stdout, ub[:,:nb], label=[label[i] for i in notimp],
                          label2=occ_label, start=1)

    mob = numpy.dot(corth[:,notimp], ub[:,:nb4cas])
    idxenv = numpy.argsort(-occb[nb4cas:]) + nb4cas
    mo_env = numpy.dot(corth[:,notimp], ub[:,idxenv])

    mocore = mo_env[:,:ncore]
    mocas = numpy.hstack((numpy.dot(corth[:,baslst],ui), mob))
    movir = mo_env[:,ncore:]

    if canonicalize or freeze_imp:
        if mf.mo_energy is None or mf.mo_coeff is None:
            fock = mf.get_hcore()
        else:
            if isinstance(mf.mo_coeff, numpy.ndarray) and mf.mo_coeff.ndim == 2:
                sc = numpy.dot(s, mf.mo_coeff)
                fock = numpy.dot(sc*mf.mo_energy, sc.T)
            else:
                sc = numpy.dot(s, mf.mo_coeff[0])
                fock = numpy.dot(sc*mf.mo_energy[0], sc.T)

        def trans(c):
            if c.shape[1] == 0:
                return c
            else:
                f1 = reduce(numpy.dot, (c.T, fock, c))
                e, u = scipy.linalg.eigh(f1)
                log.debug1('Fock eig %s', e)
                return symmetrize(mol, e, numpy.dot(c, u), s, log)

        if freeze_imp:
            log.debug('Semi-canonicalization for freeze_imp=True')
            mo = numpy.hstack([trans(mocore), trans(mocas[:,:nimp]),
                               trans(mocas[:,nimp:]), trans(movir)])
        else:
            mo = numpy.hstack([trans(mocore), trans(mocas), trans(movir)])
    else:
        mo = numpy.hstack((mocore, mocas, movir))

    return ncas, nelecas, mo

def analyze(casscf, mo_coeff=None, ci=None, verbose=logger.INFO):
    from pyscf.tools import dump_mat
    from pyscf.mcscf import addons

    if mo_coeff is None:
        mo_coeff = casscf.mo_coeff
    if ci is None:
        ci = casscf.ci
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(casscf.stdout, verbose)
    nelecas = casscf.nelecas
    ncas = casscf.ncas
    ncore = casscf.ncore
    nocc = ncore + ncas
    label = casscf.mol.spheric_labels(True)

    if hasattr(casscf.fcisolver, "make_rdm1s"):
        casdm1a, casdm1b = casscf.fcisolver.make_rdm1s(ci, ncas, nelecas)
        casdm1 = casdm1a + casdm1b
        mocore = mo_coeff[:, :ncore]
        mocas = mo_coeff[:, ncore:nocc]
        dm1b = numpy.dot(mocore, mocore.T)
        dm1a = dm1b + reduce(numpy.dot, (mocas, casdm1a, mocas.T))
        dm1b += reduce(numpy.dot, (mocas, casdm1b, mocas.T))
        dm1 = dm1a + dm1b
        if log.verbose >= logger.DEBUG1:
            log.info("alpha density matrix (on AO)")
            dump_mat.dump_tri(log.stdout, dm1a, label)
            log.info("beta density matrix (on AO)")
            dump_mat.dump_tri(log.stdout, dm1b, label)
    else:
        casdm1 = casscf.fcisolver.make_rdm1(ci, ncas, nelecas)
        mocore = mo_coeff[:, :ncore]
        mocas = mo_coeff[:, ncore:nocc]
        dm1a = numpy.dot(mocore, mocore.T) * 2 + reduce(numpy.dot, (mocas, casdm1, mocas.T))
        dm1b = None
        dm1 = dm1a

    if log.verbose >= logger.INFO:
        # note the last two args of ._eig for mc1step_symm
        occ, ucas = casscf._eig(-casdm1, ncore, nocc)
        log.info("Natural occ %s", str(-occ))
        for i, k in enumerate(numpy.argmax(abs(ucas), axis=0)):
            if ucas[k, i] < 0:
                ucas[:, i] *= -1
        mo_cas = numpy.dot(mo_coeff[:, ncore:nocc], ucas)
        log.info("Natural orbital in CAS space")
        dump_mat.dump_rec(log.stdout, mo_cas, label, start=1)

        if casscf._scf.mo_coeff is not None:
            s = reduce(numpy.dot, (casscf.mo_coeff.T, casscf._scf.get_ovlp(), casscf._scf.mo_coeff))
            idx = numpy.argwhere(abs(s) > 0.4)
            for i, j in idx:
                log.info("<mo-mcscf|mo-hf> %d  %d  %12.8f", i + 1, j + 1, s[i, j])

        if ci is not None:
            log.info("** Largest CI components **")
            log.info(" string alpha, string beta, CI coefficients")
            for c, ia, ib in fci.addons.large_ci(ci, casscf.ncas, casscf.nelecas):
                log.info("  %9s    %9s    %.12f", ia, ib, c)

        s = casscf._scf.get_ovlp()
        # casscf._scf.mulliken_pop(casscf.mol, dm1, s, verbose=log)
        casscf._scf.mulliken_pop_meta_lowdin_ao(casscf.mol, dm1, verbose=log)
    return dm1a, dm1b

    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)

def analyze(mf, verbose=logger.DEBUG):
    from pyscf.tools import dump_mat

    mo_energy = mf.mo_energy
    mo_occ = mf.mo_occ
    mo_coeff = mf.mo_coeff
    mol = mf.mol
    log = pyscf.lib.logger.Logger(mf.stdout, verbose)
    nirrep = len(mol.irrep_id)
    ovlp_ao = mf.get_ovlp()
    orbsyma = pyscf.symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb, mo_coeff[0], s=ovlp_ao)
    orbsymb = pyscf.symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb, mo_coeff[1], s=ovlp_ao)
    orbsyma = numpy.array(orbsyma)
    orbsymb = numpy.array(orbsymb)
    tot_sym = 0
    noccsa = [sum(orbsyma[mo_occ[0] > 0] == ir) for ir in mol.irrep_id]
    noccsb = [sum(orbsymb[mo_occ[1] > 0] == ir) for ir in mol.irrep_id]
    for ir in range(nirrep):
        if (noccsa[ir] + noccsb[ir]) % 2:
            tot_sym ^= mol.irrep_id[ir]
    if mol.groupname in ("Dooh", "Coov", "SO3"):
        log.info("TODO: total symmetry for %s", mol.groupname)
    else:
        log.info("total symmetry = %s", pyscf.symm.irrep_id2name(mol.groupname, tot_sym))
    log.info("alpha occupancy for each irrep:  " + (" %4s" * nirrep), *mol.irrep_name)
    log.info("                                 " + (" %4d" * nirrep), *noccsa)
    log.info("beta  occupancy for each irrep:  " + (" %4s" * nirrep), *mol.irrep_name)
    log.info("                                 " + (" %4d" * nirrep), *noccsb)

    ss, s = mf.spin_square((mo_coeff[0][:, mo_occ[0] > 0], mo_coeff[1][:, mo_occ[1] > 0]), ovlp_ao)
    log.info("multiplicity <S^2> = %.8g  2S+1 = %.8g", ss, s)

    if verbose >= logger.INFO:
        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(orbsyma):
            if j in irorbcnt:
                irorbcnt[j] += 1
            else:
                irorbcnt[j] = 1
            log.info(
                "alpha MO #%d (%s #%d), energy= %.15g occ= %g",
                k + 1,
                irname_full[j],
                irorbcnt[j],
                mo_energy[0][k],
                mo_occ[0][k],
            )
        irorbcnt = {}
        for k, j in enumerate(orbsymb):
            if j in irorbcnt:
                irorbcnt[j] += 1
            else:
                irorbcnt[j] = 1
            log.info(
                "beta  MO #%d (%s #%d), energy= %.15g occ= %g",
                k + 1,
                irname_full[j],
                irorbcnt[j],
                mo_energy[1][k],
                mo_occ[1][k],
            )

    if mf.verbose >= logger.DEBUG:
        label = mol.spheric_labels(True)
        molabel = []
        irorbcnt = {}
        for k, j in enumerate(orbsyma):
            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(" ** alpha MO coefficients **")
        dump_mat.dump_rec(mol.stdout, mo_coeff[0], label, molabel, start=1)

        molabel = []
        irorbcnt = {}
        for k, j in enumerate(orbsymb):
            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(" ** beta MO coefficients **")
        dump_mat.dump_rec(mol.stdout, mo_coeff[1], label, molabel, start=1)

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