def test_orth_ao(self):
     c0 = orth.pre_orth_ao(mol, method='scf')
     self.assertAlmostEqual(numpy.linalg.norm(c0), 5.7742626195362039, 9)
     self.assertAlmostEqual(abs(c0).sum(), 33.461490657433551, 8)
     c = orth.orth_ao(mol, 'lowdin', c0)
     self.assertAlmostEqual(numpy.linalg.norm(c), 8.9823854843222257, 9)
     self.assertAlmostEqual(abs(c).sum(), 94.933979307106767, 8)
     c = orth.orth_ao(mol, 'meta_lowdin', c0)
     self.assertAlmostEqual(numpy.linalg.norm(c), 8.9823854843222257, 9)
     self.assertAlmostEqual(abs(c).sum(), 93.029386338534394, 8)

 def test_orth_ao(self):
     c0 = orth.pre_orth_ao(mol, method='scf')
     self.assertAlmostEqual(numpy.linalg.norm(c0), 7.2617698799320358, 9)
     self.assertAlmostEqual(abs(c0).sum(), 40.116080631662804, 8)
     c = orth.orth_ao(mol, 'lowdin', c0)
     self.assertAlmostEqual(numpy.linalg.norm(c), 10.967144073462256, 9)
     self.assertAlmostEqual(abs(c).sum(), 112.23459140302003, 8)
     c = orth.orth_ao(mol, 'meta_lowdin', c0)
     self.assertAlmostEqual(numpy.linalg.norm(c), 10.967144073462256, 9)
     self.assertAlmostEqual(abs(c).sum(), 111.61017124719302, 8)

def mulliken_meta(cell, dm_ao_kpts, verbose=logger.DEBUG,
                  pre_orth_method=PRE_ORTH_METHOD, s=None):
    '''A modified Mulliken population analysis, based on meta-Lowdin AOs.

    Note this function only computes the Mulliken population for the gamma
    point density matrix.
    '''
    from pyscf.lo import orth
    if s is None:
        s = khf.get_ovlp(cell)
    log = logger.new_logger(cell, verbose)
    log.note('Analyze output for *gamma point*.')
    log.info('    To include the contributions from k-points, transform to a '
             'supercell then run the population analysis on the supercell\n'
             '        from pyscf.pbc.tools import k2gamma\n'
             '        k2gamma.k2gamma(mf).mulliken_meta()')
    log.note("KUHF mulliken_meta")
    dm_ao_gamma = dm_ao_kpts[:,0,:,:].real
    s_gamma = s[0,:,:].real
    c = orth.restore_ao_character(cell, pre_orth_method)
    orth_coeff = orth.orth_ao(cell, 'meta_lowdin', pre_orth_ao=c, s=s_gamma)
    c_inv = np.dot(orth_coeff.T, s_gamma)
    dm_a = reduce(np.dot, (c_inv, dm_ao_gamma[0], c_inv.T.conj()))
    dm_b = reduce(np.dot, (c_inv, dm_ao_gamma[1], c_inv.T.conj()))

    log.note(' ** Mulliken pop alpha/beta on meta-lowdin orthogonal AOs **')
    return mol_uhf.mulliken_pop(cell, (dm_a,dm_b), np.eye(orth_coeff.shape[0]), log)

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 test_ghost_atm_meta_lowdin(self):
     mol = gto.Mole()
     mol.atom = [["O" , (0. , 0.     , 0.)],
                 ['ghost'   , (0. , -0.757, 0.587)],
                 [1   , (0. , 0.757 , 0.587)] ]
     mol.spin = 1
     mol.basis = {'O':'ccpvdz', 'H':'ccpvdz',
                  'GHOST': gto.basis.load('631g','H')}
     mol.build()
     c = orth.orth_ao(mol, method='meta_lowdin')
     self.assertAlmostEqual(numpy.linalg.norm(c), 7.9067188905237256, 9)

def atomic_init_guess(mol, mo_coeff):
    s = mol.intor_symmetric('int1e_ovlp')
    c = orth.orth_ao(mol, s=s)
    mo = reduce(numpy.dot, (c.conj().T, s, mo_coeff))
    nmo = mo_coeff.shape[1]
# Find the AOs which have largest overlap to MOs
    idx = numpy.argsort(numpy.einsum('pi,pi->p', mo.conj(), mo))
    nmo = mo.shape[1]
    idx = idx[-nmo:]
    u, w, vh = numpy.linalg.svd(mo[idx])
    return lib.dot(vh, u.conj().T)

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 atomic_init_guess(mol, mo_coeff):
    s = mol.intor_symmetric('int1e_ovlp')
    c = orth.orth_ao(mol, s=s)
    mo = reduce(numpy.dot, (c.conj().T, s, mo_coeff))
# Find the AOs which have largest overlap to MOs
    idx = numpy.argsort(numpy.einsum('pi,pi->p', mo.conj(), mo))
    nmo = mo.shape[1]
    idx = sorted(idx[-nmo:])

    # Rotate mo_coeff, make it as close as possible to AOs
    u, w, vh = numpy.linalg.svd(mo[idx])
    return lib.dot(u, vh).conj().T

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 mulliken_meta(mol, dm_ao, verbose=logger.DEBUG,
                  pre_orth_method=PRE_ORTH_METHOD, s=None):
    '''Mulliken population analysis, based on meta-Lowdin AOs.
    '''
    from pyscf.lo import orth
    if s is None: s = hf.get_ovlp(mol)
    log = logger.new_logger(mol, verbose)
    if isinstance(dm_ao, numpy.ndarray) and dm_ao.ndim == 2:
        dm_ao = numpy.array((dm_ao*.5, dm_ao*.5))
    c = orth.restore_ao_character(mol, pre_orth_method)
    orth_coeff = orth.orth_ao(mol, 'meta_lowdin', pre_orth_ao=c, s=s)
    c_inv = numpy.dot(orth_coeff.T, s)
    dm_a = reduce(numpy.dot, (c_inv, dm_ao[0], c_inv.T.conj()))
    dm_b = reduce(numpy.dot, (c_inv, dm_ao[1], c_inv.T.conj()))

    log.note(' ** Mulliken pop alpha/beta on meta-lowdin orthogonal AOs **')
    return mulliken_pop(mol, (dm_a,dm_b), numpy.eye(orth_coeff.shape[0]), log)

def mulliken_pop_meta_lowdin_ao(mol, dm_ao, verbose=logger.DEBUG,
                                pre_orth_method='ANO'):
    '''Mulliken population analysis, based on meta-Lowdin AOs.
    '''
    from pyscf.lo import orth
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(mol.stdout, verbose)
    c = orth.pre_orth_ao(mol, pre_orth_method)
    orth_coeff = orth.orth_ao(mol, 'meta_lowdin', pre_orth_ao=c)
    c_inv = numpy.linalg.inv(orth_coeff)
    dm_a = reduce(numpy.dot, (c_inv, dm_ao[0], c_inv.T.conj()))
    dm_b = reduce(numpy.dot, (c_inv, dm_ao[1], c_inv.T.conj()))

    log.info(' ** Mulliken pop alpha/beta on meta-lowdin orthogonal AOs **')
    return mulliken_pop(mol, (dm_a,dm_b), numpy.eye(orth_coeff.shape[0]), 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 mulliken_meta(mol, dm, verbose=logger.DEBUG, pre_orth_method='ANO',
                  s=None):
    '''Mulliken population analysis, based on meta-Lowdin AOs.
    In the meta-lowdin, the AOs are grouped in three sets: core, valence and
    Rydberg, the orthogonalization are carreid out within each subsets.

    Args:
        mol : an instance of :class:`Mole`

        dm : ndarray or 2-item list of ndarray
            Density matrix.  ROHF dm is a 2-item list of 2D array

    Kwargs:
        verbose : int or instance of :class:`lib.logger.Logger`

        pre_orth_method : str
            Pre-orthogonalization, which localized GTOs for each atom.
            To obtain the occupied and unoccupied atomic shells, there are
            three methods

            | 'ano'   : Project GTOs to ANO basis
            | 'minao' : Project GTOs to MINAO basis
            | 'scf'   : Fraction-averaged RHF

    '''
    from pyscf.lo import orth
    if s is None:
        s = get_ovlp(mol)
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(mol.stdout, verbose)
    c = orth.pre_orth_ao(mol, pre_orth_method)
    orth_coeff = orth.orth_ao(mol, 'meta_lowdin', pre_orth_ao=c, s=s)
    c_inv = numpy.dot(orth_coeff.T, s)
    if isinstance(dm, numpy.ndarray) and dm.ndim == 2:
        dm = reduce(numpy.dot, (c_inv, dm, c_inv.T.conj()))
    else:  # ROHF
        dm = reduce(numpy.dot, (c_inv, dm[0]+dm[1], c_inv.T.conj()))

    log.info(' ** Mulliken pop on meta-lowdin orthogonal AOs  **')
    return mulliken_pop(mol, dm, numpy.eye(orth_coeff.shape[0]), log)

def atomic_pops(mol, mo_coeff, method='meta_lowdin'):
    '''kwarg method can be one of mulliken, lowdin, meta_lowdin
    '''
    s = mol.intor_symmetric('cint1e_ovlp_sph')
    nmo = mo_coeff.shape[1]
    proj = numpy.empty((mol.natm,nmo,nmo))

    if method.lower() == 'mulliken':
        for i, (b0, b1, p0, p1) in enumerate(mol.offset_nr_by_atom()):
            csc = reduce(numpy.dot, (mo_coeff[p0:p1].T, s[p0:p1], mo_coeff))
            proj[i] = (csc + csc.T) * .5

    elif method.lower() in ('lowdin', 'meta_lowdin'):
        csc = reduce(numpy.dot, (mo_coeff.T, s, orth.orth_ao(mol, method, s=s)))
        for i, (b0, b1, p0, p1) in enumerate(mol.offset_nr_by_atom()):
            proj[i] = numpy.dot(csc[:,p0:p1], csc[:,p0:p1].T)
    else:
        raise KeyError('method = %s' % method)

    return proj

def mulliken_meta(mol, dm_ao, verbose=logger.DEBUG, pre_orth_method="ANO", s=None):
    """Mulliken population analysis, based on meta-Lowdin AOs.
    """
    from pyscf.lo import orth

    if s is None:
        s = hf.get_ovlp(mol)
    if isinstance(verbose, logger.Logger):
        log = verbose
    else:
        log = logger.Logger(mol.stdout, verbose)
    if isinstance(dm_ao, numpy.ndarray) and dm_ao.ndim == 2:
        dm_ao = numpy.array((dm_ao * 0.5, dm_ao * 0.5))
    c = orth.pre_orth_ao(mol, pre_orth_method)
    orth_coeff = orth.orth_ao(mol, "meta_lowdin", pre_orth_ao=c, s=s)
    c_inv = numpy.dot(orth_coeff.T, s)
    dm_a = reduce(numpy.dot, (c_inv, dm_ao[0], c_inv.T.conj()))
    dm_b = reduce(numpy.dot, (c_inv, dm_ao[1], c_inv.T.conj()))

    log.info(" ** Mulliken pop alpha/beta on meta-lowdin orthogonal AOs **")
    return mulliken_pop(mol, (dm_a, dm_b), numpy.eye(orth_coeff.shape[0]), log)

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(casscf, mo_coeff=None, ci=None, verbose=None,
            large_ci_tol=LARGE_CI_TOL, with_meta_lowdin=WITH_META_LOWDIN,
            **kwargs):
    from pyscf.lo import orth
    from pyscf.tools import dump_mat
    from pyscf.mcscf import addons
    log = logger.new_logger(casscf, verbose)

    if mo_coeff is None: mo_coeff = casscf.mo_coeff
    if ci is None: ci = casscf.ci
    nelecas = casscf.nelecas
    ncas = casscf.ncas
    ncore = casscf.ncore
    nocc = ncore + ncas
    mocore = mo_coeff[:,:ncore]
    mocas = mo_coeff[:,ncore:nocc]

    label = casscf.mol.ao_labels()
    if (isinstance(ci, (list, tuple)) and
        not isinstance(casscf.fcisolver, addons.StateAverageFCISolver)):
        log.warn('Mulitple states found in CASCI/CASSCF solver. Density '
                 'matrix of first state is generated in .analyze() function.')
        civec = ci[0]
    else:
        civec = ci
    if getattr(casscf.fcisolver, 'make_rdm1s', None):
        casdm1a, casdm1b = casscf.fcisolver.make_rdm1s(civec, ncas, nelecas)
        casdm1 = casdm1a + casdm1b
        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.DEBUG2:
            log.info('alpha density matrix (on AO)')
            dump_mat.dump_tri(log.stdout, dm1a, label, **kwargs)
            log.info('beta density matrix (on AO)')
            dump_mat.dump_tri(log.stdout, dm1b, label, **kwargs)
    else:
        casdm1 = casscf.fcisolver.make_rdm1(civec, ncas, nelecas)
        dm1a =(numpy.dot(mocore, mocore.T) * 2
             + reduce(numpy.dot, (mocas, casdm1, mocas.T)))
        dm1b = None
        dm1 = dm1a

    if log.verbose >= logger.INFO:
        ovlp_ao = casscf._scf.get_ovlp()
        # note the last two args of ._eig for mc1step_symm
        occ, ucas = casscf._eig(-casdm1, ncore, nocc)
        log.info('Natural occ %s', str(-occ))
        mocas = numpy.dot(mocas, ucas)
        if with_meta_lowdin:
            log.info('Natural orbital (expansion on meta-Lowdin AOs) in CAS space')
            orth_coeff = orth.orth_ao(casscf.mol, 'meta_lowdin', s=ovlp_ao)
            mocas = reduce(numpy.dot, (orth_coeff.T, ovlp_ao, mocas))
        else:
            log.info('Natural orbital (expansion on AOs) in CAS space')
        dump_mat.dump_rec(log.stdout, mocas, label, start=1, **kwargs)
        if log.verbose >= logger.DEBUG2:
            if not casscf.natorb:
                log.debug2('NOTE: mc.mo_coeff in active space is different to '
                           'the natural orbital coefficients printed in above.')
            if with_meta_lowdin:
                c = reduce(numpy.dot, (orth_coeff.T, ovlp_ao, mo_coeff))
                log.debug2('MCSCF orbital (expansion on meta-Lowdin AOs)')
            else:
                c = mo_coeff
                log.debug2('MCSCF orbital (expansion on AOs)')
            dump_mat.dump_rec(log.stdout, c, label, start=1, **kwargs)

        if casscf._scf.mo_coeff is not None:
            addons.map2hf(casscf, casscf._scf.mo_coeff)

        if getattr(casscf.fcisolver, 'large_ci', None) and ci is not None:
            log.info('** Largest CI components **')
            if isinstance(ci, (tuple, list)):
                for i, civec in enumerate(ci):
                    res = casscf.fcisolver.large_ci(civec, casscf.ncas, casscf.nelecas,
                                                    large_ci_tol, return_strs=False)
                    log.info('  [alpha occ-orbitals] [beta occ-orbitals]  state %-3d CI coefficient', i)
                    for c,ia,ib in res:
                        log.info('  %-20s %-30s %.12f', ia, ib, c)
            else:
                log.info('  [alpha occ-orbitals] [beta occ-orbitals]            CI coefficient')
                res = casscf.fcisolver.large_ci(ci, casscf.ncas, casscf.nelecas,
                                                large_ci_tol, return_strs=False)
                for c,ia,ib in res:
                    log.info('  %-20s %-30s %.12f', ia, ib, c)

        if with_meta_lowdin:
            casscf._scf.mulliken_meta(casscf.mol, dm1, s=ovlp_ao, verbose=log)
        else:
            casscf._scf.mulliken_pop(casscf.mol, dm1, s=ovlp_ao, 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.lo import orth
    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 = symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb,
                                  mo_coeff[0], s=ovlp_ao)
    orbsymb = 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 i, ir in enumerate(mol.irrep_id):
        if (noccsa[i]+noccsb[i]) % 2:
            tot_sym ^= ir
    if mol.groupname in ('Dooh', 'Coov', 'SO3'):
        log.note('TODO: total symmetry for %s', mol.groupname)
    else:
        log.note('total symmetry = %s',
                 symm.irrep_id2name(mol.groupname, tot_sym))
    log.note('alpha occupancy for each irrep:  '+(' %4s'*nirrep),
             *mol.irrep_name)
    log.note('                                 '+(' %4d'*nirrep),
             *noccsa)
    log.note('beta  occupancy for each irrep:  '+(' %4s'*nirrep),
             *mol.irrep_name)
    log.note('                                 '+(' %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.note('multiplicity <S^2> = %.8g  2S+1 = %.8g', ss, s)

    if verbose >= logger.NOTE:
        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(orbsyma):
            if j in irorbcnt:
                irorbcnt[j] += 1
            else:
                irorbcnt[j] = 1
            log.note('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.note('beta  MO #%d (%s #%d), energy= %.15g occ= %g',
                     k+1, irname_full[j], irorbcnt[j], mo_energy[1][k], mo_occ[1][k])

    ovlp_ao = mf.get_ovlp()
    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 (expansion on meta-Lowdin AOs) **')
        orth_coeff = orth.orth_ao(mol, 'meta_lowdin', s=ovlp_ao)
        c_inv = numpy.dot(orth_coeff.T, ovlp_ao)
        dump_mat.dump_rec(mol.stdout, c_inv.dot(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 (expansion on meta-Lowdin AOs) **')
        dump_mat.dump_rec(mol.stdout, c_inv.dot(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)