def _get_basis_type(mol):
    def classify(mol_basis):
        basis_type = 'other'
        if isinstance(mol_basis, (str, unicode)):
            mol_basis = gto.basis._format_basis_name(mol_basis)
            if mol_basis[:6] == 'def2tz':
                basis_type = 'def2-TZ'
            elif mol_basis[:6] == 'def2sv':
                basis_type = 'sv'
            elif mol_basis[:5] == '631g*':
                basis_type = '6-31gd'
            elif mol_basis[:4] == '631g' and 'd' in mol_basis:
                basis_type = '6-31gd'
        return basis_type

    if isinstance(mol.basis, dict):
        basis_types = [classify(b) for b in mol.basis.values()]
        basis_type = 'other'
        for bt in basis_types:
            if bt != 'other':
                basis_type = bt
                break
        if (len(basis_types) > 1 and
            all(b == basis_type for b in basis_types)):
            logger.warn('Mutliple types of basis found in mol.basis. '
                        'Type %s is applied\n' % basis_type)
    else:
        basis_type = classify(mol.basis)
    return basis_type

    def _common_init_(self, mycc, mo_coeff=None):
        if mo_coeff is None:
            mo_coeff = mycc.mo_coeff
        mo_idx = ccsd.get_frozen_mask(mycc)
        if getattr(mo_coeff, 'orbspin', None) is not None:
            self.orbspin = mo_coeff.orbspin[mo_idx]
            mo_coeff = lib.tag_array(mo_coeff[:,mo_idx], orbspin=self.orbspin)
            self.mo_coeff = mo_coeff
        else:
            orbspin = scf.ghf.guess_orbspin(mo_coeff)
            self.mo_coeff = mo_coeff = mo_coeff[:,mo_idx]
            if not np.any(orbspin == -1):
                self.orbspin = orbspin[mo_idx]
                self.mo_coeff = lib.tag_array(mo_coeff, orbspin=self.orbspin)

        # Note: Recomputed fock matrix since SCF may not be fully converged.
        dm = mycc._scf.make_rdm1(mycc.mo_coeff, mycc.mo_occ)
        fockao = mycc._scf.get_fock(dm=dm)
        self.fock = reduce(np.dot, (mo_coeff.conj().T, fockao, mo_coeff))
        self.nocc = mycc.nocc

        mo_e = self.mo_energy = self.fock.diagonal().real
        gap = abs(mo_e[:self.nocc,None] - mo_e[None,self.nocc:]).min()
        if gap < 1e-5:
            logger.warn(mycc, 'HOMO-LUMO gap %s too small for GCCSD', gap)
        return self

def _nao_sub(mol, pre_occ, pre_nao, s=None):
    if s is None:
        s = mol.intor_symmetric('cint1e_ovlp_sph')
    core_lst, val_lst, rydbg_lst = _core_val_ryd_list(mol)
    nbf = mol.nao_nr()
    cnao = numpy.empty((nbf,nbf))

    if core_lst:
        c = pre_nao[:,core_lst]
        s1 = reduce(numpy.dot, (c.T, s, c))
        cnao[:,core_lst] = c1 = numpy.dot(c, orth.lowdin(s1))
        c = pre_nao[:,val_lst]
        c -= reduce(numpy.dot, (c1, c1.T, s, c))
    else:
        c = pre_nao[:,val_lst]

    s1 = reduce(numpy.dot, (c.T, s, c))
    wt = pre_occ[val_lst]
    cnao[:,val_lst] = numpy.dot(c, orth.weight_orth(s1, wt))

    if rydbg_lst:
        cvlst = core_lst + val_lst
        c1 = cnao[:,cvlst]
        c = pre_nao[:,rydbg_lst]
        c -= reduce(numpy.dot, (c1, c1.T, s, c))
        s1 = reduce(numpy.dot, (c.T, s, c))
        cnao[:,rydbg_lst] = numpy.dot(c, orth.lowdin(s1))
    snorm = numpy.linalg.norm(reduce(numpy.dot, (cnao.T, s, cnao)) - numpy.eye(nbf))
    if snorm > 1e-9:
        logger.warn(mol, 'Weak orthogonality for localized orbitals %s', snorm)
    return cnao

def label_orb_symm(mol, irrep_name, symm_orb, mo, s=None, check=True):
    '''Label the symmetry of given orbitals

    irrep_name can be either the symbol or the ID of the irreducible
    representation.  If the ID is provided, it returns the numeric code
    associated with XOR operator, see :py:meth:`symm.param.IRREP_ID_TABLE`

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

        irrep_name : list of str or int
            A list of irrep ID or name,  it can be either mol.irrep_id or
            mol.irrep_name.  It can affect the return "label".
        symm_orb : list of 2d array
            the symmetry adapted basis
        mo : 2d array
            the orbitals to label

    Returns:
        list of symbols or integers to represent the irreps for the given
        orbitals

    Examples:

    >>> from pyscf import gto, scf, symm
    >>> mol = gto.M(atom='H 0 0 0; H 0 0 1', basis='ccpvdz',verbose=0, symmetry=1)
    >>> mf = scf.RHF(mol)
    >>> mf.kernel()
    >>> symm.label_orb_symm(mol, mol.irrep_name, mol.symm_orb, mf.mo_coeff)
    ['Ag', 'B1u', 'Ag', 'B1u', 'B2u', 'B3u', 'Ag', 'B2g', 'B3g', 'B1u']
    >>> symm.label_orb_symm(mol, mol.irrep_id, mol.symm_orb, mf.mo_coeff)
    [0, 5, 0, 5, 6, 7, 0, 2, 3, 5]
    '''
    nmo = mo.shape[1]
    if s is None:
        s = mol.intor_symmetric('cint1e_ovlp_sph')
    mo_s = numpy.dot(mo.T, s)
    norm = numpy.empty((len(irrep_name), nmo))
    for i,ir in enumerate(irrep_name):
        moso = numpy.dot(mo_s, symm_orb[i])
        norm[i] = numpy.einsum('ij,ij->i', moso, moso)
    iridx = numpy.argmax(norm, axis=0)
    orbsym = [irrep_name[i] for i in iridx]
    logger.debug(mol, 'irreps of each MO %s', str(orbsym))
    if check:
        norm[iridx,numpy.arange(nmo)] = 0
        orbidx = numpy.where(norm > THRESHOLD)
        if orbidx[1].size > 0:
            idx = numpy.where(norm > THRESHOLD*1e2)
            if idx[1].size > 0:
                logger.error(mol, 'orbitals %s not symmetrized, norm = %s',
                             idx[1], norm[idx])
                raise ValueError('orbitals %s not symmetrized' % idx[1])
            else:
                logger.warn(mol, 'orbitals %s not strictly symmetrized.',
                            orbidx[1])
                logger.warn(mol, 'They can be symmetrized with '
                            'pyscf.symm.symmetrize_orb function.')
                logger.debug(mol, 'norm = %s', norm[orbidx])
    return orbsym

 def get_init_guess(self, mol=None, key='minao'):
     if mol is None:
         mol = self.mol
     if callable(key):
         dm = key(mol)
     elif key.lower() == '1e':
         dm = self.init_guess_by_1e(mol)
     elif getattr(mol, 'natm', 0) == 0:
         logger.info(self, 'No atom found in mol. Use 1e initial guess')
         dm = self.init_guess_by_1e(mol)
     elif key.lower() == 'atom':
         dm = self.init_guess_by_atom(mol)
     elif key.lower() == 'chkfile':
         try:
             dm = self.init_guess_by_chkfile()
         except (IOError, KeyError):
             logger.warn(self, 'Fail in reading %s. Use MINAO initial guess',
                         self.chkfile)
             dm = self.init_guess_by_minao(mol)
     else:
         dm = self.init_guess_by_minao(mol)
     if self.verbose >= logger.DEBUG1:
         logger.debug1(self, 'Nelec from initial guess = %g',
                       (dm*self.get_ovlp()).sum().real)
     return dm

    def get_occ(self, mo_energy=None, mo_coeff=None):
        '''Label the occupancies for each orbital

        Kwargs:
            mo_energy : 1D ndarray
                Obital energies

            mo_coeff : 2D ndarray
                Obital coefficients

        Examples:

        >>> from pyscf import gto, scf
        >>> mol = gto.M(atom='H 0 0 0; F 0 0 1.1')
        >>> mf = scf.hf.SCF(mol)
        >>> mf.get_occ(numpy.arange(mol.nao_nr()))
        array([2, 2, 2, 2, 2, 0])
        '''
        if mo_energy is None: mo_energy = self.mo_energy
        mo_occ = numpy.zeros_like(mo_energy)
        nocc = self.mol.nelectron // 2
        mo_occ[:nocc] = 2
        if nocc < mo_occ.size:
            logger.info(self, 'HOMO = %.12g, LUMO = %.12g,',
                        mo_energy[nocc-1], mo_energy[nocc])
            if mo_energy[nocc-1]+1e-3 > mo_energy[nocc]:
                logger.warn(self, '!! HOMO %.12g == LUMO %.12g',
                            mo_energy[nocc-1], mo_energy[nocc])
        else:
            logger.info(self, 'HOMO = %.12g,', mo_energy[nocc-1])
        if self.verbose >= logger.DEBUG:
            numpy.set_printoptions(threshold=len(mo_energy))
            logger.debug(self, '  mo_energy = %s', mo_energy)
            numpy.set_printoptions()
        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 makov_payne_correction(mf):
    '''Makov-Payne correction (Phys. Rev. B, 51, 4014)
    '''
    cell = mf.cell
    logger.note(mf, 'Makov-Payne correction for charged 3D PBC systems')
    # PRB 51 (1995), 4014
    # PRB 77 (2008), 115139
    if cell.dimension != 3:
        logger.warn(mf, 'Correction for low-dimension PBC systems'
                    'is not available.')
        return 0

    de_mono, de_dip, de_quad, de = _dip_correction(mf)

    if mf.verbose >= logger.NOTE:
        write = mf.stdout.write
        write('Corrections (AU)\n')
        write('       Monopole      Dipole          Quadrupole    total\n')
        write('SC   %12.8f   %12.8f   %12.8f   %12.8f\n' %
              (de_mono[0], de_dip   , de_quad   , de[0]))
        write('BCC  %12.8f   %12.8f   %12.8f   %12.8f\n' %
              (de_mono[1], de_dip   , de_quad   , de[1]))
        write('FCC  %12.8f   %12.8f   %12.8f   %12.8f\n' %
              (de_mono[2], de_dip   , de_quad   , de[2]))
    return de

    def get_occ(mo_energy_kpts=None, mo_coeff=None):
        if mo_energy_kpts is None: mo_energy_kpts = mf.mo_energy

        if nelec is None:
            cell_nelec = mf.cell.nelec
        else:
            cell_nelec = nelec

        homo=[-1e8,-1e8]
        lumo=[1e8,1e8]
        mo_occ_kpts = [[], []]
        for s in [0,1]:
            for k, mo_energy in enumerate(mo_energy_kpts[s]):
                e_idx = numpy.argsort(mo_energy)
                e_sort = mo_energy[e_idx]
                n = cell_nelec[s]
                mo_occ = numpy.zeros_like(mo_energy)
                mo_occ[e_idx[:n]] = 1
                homo[s] = max(homo[s], e_sort[n-1])
                lumo[s] = min(lumo[s], e_sort[n])
                mo_occ_kpts[s].append(mo_occ)

        for nm,s in zip(['alpha','beta'],[0,1]):
            logger.info(mf, nm+' HOMO = %.12g  LUMO = %.12g', homo[s], lumo[s])
            if homo[s] > lumo[s]:
                logger.warn(mf, "WARNING! HOMO is greater than LUMO! "
                            "This may lead to incorrect canonical occupation.")

        return mo_occ_kpts

def get_occ(mf, mo_energy=None, mo_coeff=None):
    '''Label the occupancies for each orbital.
    NOTE the occupancies are not assigned based on the orbital energy ordering.
    The first N orbitals are assigned to be occupied orbitals.

    Examples:

    >>> mol = gto.M(atom='H 0 0 0; O 0 0 1.1', spin=1)
    >>> mf = scf.hf.SCF(mol)
    >>> energy = numpy.array([-10., -1., 1, -2., 0, -3])
    >>> mf.get_occ(energy)
    array([2, 2, 2, 2, 1, 0])
    '''

    if mo_energy is None: mo_energy = mf.mo_energy
    if getattr(mo_energy, 'mo_ea', None) is not None:
        mo_ea = mo_energy.mo_ea
        mo_eb = mo_energy.mo_eb
    else:
        mo_ea = mo_eb = mo_energy
    nmo = mo_ea.size
    mo_occ = numpy.zeros(nmo)
    if getattr(mf, 'nelec', None) is None:
        nelec = mf.mol.nelec
    else:
        nelec = mf.nelec
    ncore = nelec[1]
    nocc  = nelec[0]
    nopen = abs(nocc - ncore)
    mo_occ = _fill_rohf_occ(mo_energy, mo_ea, mo_eb, ncore, nopen)

    if mf.verbose >= logger.INFO and nocc < nmo and ncore > 0:
        ehomo = max(mo_energy[mo_occ> 0])
        elumo = min(mo_energy[mo_occ==0])
        if ehomo+1e-3 > elumo:
            logger.warn(mf, 'HOMO %.15g >= LUMO %.15g', ehomo, elumo)
        else:
            logger.info(mf, '  HOMO = %.15g  LUMO = %.15g', ehomo, elumo)
        if nopen > 0 and mf.verbose >= logger.DEBUG:
            core_idx = mo_occ == 2
            open_idx = mo_occ == 1
            vir_idx = mo_occ == 0
            logger.debug(mf, '                  Roothaan           | alpha              | beta')
            logger.debug(mf, '  Highest 2-occ = %18.15g | %18.15g | %18.15g',
                         max(mo_energy[core_idx]),
                         max(mo_ea[core_idx]), max(mo_eb[core_idx]))
            logger.debug(mf, '  Lowest 0-occ =  %18.15g | %18.15g | %18.15g',
                         min(mo_energy[vir_idx]),
                         min(mo_ea[vir_idx]), min(mo_eb[vir_idx]))
            for i in numpy.where(open_idx)[0]:
                logger.debug(mf, '  1-occ =         %18.15g | %18.15g | %18.15g',
                             mo_energy[i], mo_ea[i], mo_eb[i])

        if mf.verbose >= logger.DEBUG:
            numpy.set_printoptions(threshold=nmo)
            logger.debug(mf, '  Roothaan mo_energy =\n%s', mo_energy)
            logger.debug1(mf, '  alpha mo_energy =\n%s', mo_ea)
            logger.debug1(mf, '  beta  mo_energy =\n%s', mo_eb)
            numpy.set_printoptions(threshold=1000)
    return mo_occ

    def _common_init_(self, mycc, mo_coeff=None):
        if mo_coeff is None:
            mo_coeff = mycc.mo_coeff
        mo_idx = mycc.get_frozen_mask()
        self.mo_coeff = mo_coeff = \
                (mo_coeff[0][:,mo_idx[0]], mo_coeff[1][:,mo_idx[1]])
# Note: Recomputed fock matrix since SCF may not be fully converged.
        dm = mycc._scf.make_rdm1(mycc.mo_coeff, mycc.mo_occ)
        fockao = mycc._scf.get_fock(dm=dm)
        self.focka = reduce(np.dot, (mo_coeff[0].conj().T, fockao[0], mo_coeff[0]))
        self.fockb = reduce(np.dot, (mo_coeff[1].conj().T, fockao[1], mo_coeff[1]))
        self.fock = (self.focka, self.fockb)
        nocca, noccb = self.nocc = mycc.nocc
        self.mol = mycc.mol

        mo_ea = self.focka.diagonal().real
        mo_eb = self.fockb.diagonal().real
        self.mo_energy = (mo_ea, mo_eb)
        gap_a = abs(mo_ea[:nocca,None] - mo_ea[None,nocca:])
        gap_b = abs(mo_eb[:noccb,None] - mo_eb[None,noccb:])
        if gap_a.size > 0:
            gap_a = gap_a.min()
        else:
            gap_a = 1e9
        if gap_b.size > 0:
            gap_b = gap_b.min()
        else:
            gap_b = 1e9
        if gap_a < 1e-5 or gap_b < 1e-5:
            logger.warn(mycc, 'HOMO-LUMO gap (%s,%s) too small for UCCSD',
                        gap_a, gap_b)
        return self

def energy(cc, t1=None, t2=None, eris=None):
    '''UCCSD correlation energy'''
    if t1 is None: t1 = cc.t1
    if t2 is None: t2 = cc.t2
    if eris is None: eris = cc.ao2mo()

    t1a, t1b = t1
    t2aa, t2ab, t2bb = t2
    nocca, noccb, nvira, nvirb = t2ab.shape
    eris_ovov = np.asarray(eris.ovov)
    eris_OVOV = np.asarray(eris.OVOV)
    eris_ovOV = np.asarray(eris.ovOV)
    fova = eris.focka[:nocca,nocca:]
    fovb = eris.fockb[:noccb,noccb:]
    e  = np.einsum('ia,ia', fova, t1a)
    e += np.einsum('ia,ia', fovb, t1b)
    e += 0.25*np.einsum('ijab,iajb',t2aa,eris_ovov)
    e -= 0.25*np.einsum('ijab,ibja',t2aa,eris_ovov)
    e += 0.25*np.einsum('ijab,iajb',t2bb,eris_OVOV)
    e -= 0.25*np.einsum('ijab,ibja',t2bb,eris_OVOV)
    e +=      np.einsum('iJaB,iaJB',t2ab,eris_ovOV)
    e += 0.5*np.einsum('ia,jb,iajb',t1a,t1a,eris_ovov)
    e -= 0.5*np.einsum('ia,jb,ibja',t1a,t1a,eris_ovov)
    e += 0.5*np.einsum('ia,jb,iajb',t1b,t1b,eris_OVOV)
    e -= 0.5*np.einsum('ia,jb,ibja',t1b,t1b,eris_OVOV)
    e +=     np.einsum('ia,jb,iajb',t1a,t1b,eris_ovOV)
    if abs(e.imag) > 1e-4:
        logger.warn(cc, 'Non-zero imaginary part found in UCCSD energy %s', e)
    return e.real

def get_fermi(mf, mo_energy_kpts=None, mo_occ_kpts=None):
    '''A pair of Fermi level for spin-up and spin-down orbitals
    '''
    if mo_energy_kpts is None: mo_energy_kpts = mf.mo_energy
    if mo_occ_kpts is None: mo_occ_kpts = mf.mo_occ

    # mo_energy_kpts and mo_occ_kpts are k-point UHF quantities
    assert(mo_energy_kpts[0][0].ndim == 1)
    assert(mo_occ_kpts[0][0].ndim == 1)

    nocca = sum(mo_occ.sum() for mo_occ in mo_occ_kpts[0])
    noccb = sum(mo_occ.sum() for mo_occ in mo_occ_kpts[1])
    # nocc may not be perfect integer when smearing is enabled
    nocca = int(nocca.round(3))
    noccb = int(noccb.round(3))

    fermi_a = np.sort(np.hstack(mo_energy_kpts[0]))[nocca-1]
    fermi_b = np.sort(np.hstack(mo_energy_kpts[1]))[noccb-1]

    for k, mo_e in enumerate(mo_energy_kpts[0]):
        mo_occ = mo_occ_kpts[0][k]
        if mo_occ[mo_e > fermi_a].sum() > 0.5:
            logger.warn(mf, 'Alpha occupied band above Fermi level: \n'
                        'k=%d, mo_e=%s, mo_occ=%s', k, mo_e, mo_occ)
    for k, mo_e in enumerate(mo_energy_kpts[1]):
        mo_occ = mo_occ_kpts[1][k]
        if mo_occ[mo_e > fermi_b].sum() > 0.5:
            logger.warn(mf, 'Beta occupied band above Fermi level: \n'
                        'k=%d, mo_e=%s, mo_occ=%s', k, mo_e, mo_occ)
    return (fermi_a, fermi_b)

    def get_occ(self, mo_energy=None, mo_coeff=None):
        '''Label the occupancies for each orbital.
        NOTE the occupancies are not assigned based on the orbital energy ordering.
        The first N orbitals are assigned to be occupied orbitals.

        Examples:

        >>> mol = gto.M(atom='H 0 0 0; O 0 0 1.1', spin=1)
        >>> mf = scf.hf.SCF(mol)
        >>> energy = numpy.array([-10., -1., 1, -2., 0, -3])
        >>> mf.get_occ(energy)
        array([2, 2, 2, 2, 1, 0])
        '''

        if mo_energy is None: mo_energy = self.mo_energy
        mo_occ = numpy.zeros_like(mo_energy)
        ncore = self.nelec[1]
        nopen = self.nelec[0] - ncore
        nocc = ncore + nopen
        mo_occ[:ncore] = 2
        mo_occ[ncore:nocc] = 1
        if nocc < len(mo_energy):
            logger.info(self, 'HOMO = %.12g  LUMO = %.12g',
                        mo_energy[nocc-1], mo_energy[nocc])
            if mo_energy[nocc-1]+1e-3 > mo_energy[nocc]:
                logger.warn(self.mol, '!! HOMO %.12g == LUMO %.12g',
                            mo_energy[nocc-1], mo_energy[nocc])
        else:
            logger.info(self, 'HOMO = %.12g  no LUMO', mo_energy[nocc-1])
        if nopen > 0:
            for i in range(ncore, nocc):
                logger.debug(self, 'singly occupied orbital energy = %.12g',
                             mo_energy[i])
        logger.debug(self, '  mo_energy = %s', mo_energy)
        return mo_occ

    def build(self, mol=None):
        if mol is None: mol = self.mol
        if mol.symmetry:
            for irname in self.irrep_nelec:
                if irname not in mol.irrep_name:
                    logger.warn(self, 'Molecule does not have irrep %s', irname)

            nelec_fix = self.irrep_nelec.values()
            if any(isinstance(x, (tuple, list)) for x in nelec_fix):
                msg =('Number of alpha/beta electrons cannot be assigned '
                      'separately in GHF.  irrep_nelec = %s' % self.irrep_nelec)
                raise ValueError(msg)
            nelec_fix = sum(nelec_fix)
            float_irname = set(mol.irrep_name) - set(self.irrep_nelec)
            if nelec_fix > mol.nelectron:
                msg =('More electrons defined by irrep_nelec than total num electrons. '
                      'mol.nelectron = %d  irrep_nelec = %s' %
                      (mol.nelectron, self.irrep_nelec))
                raise ValueError(msg)
            else:
                logger.info(mol, 'Freeze %d electrons in irreps %s',
                            nelec_fix, self.irrep_nelec.keys())

            if len(float_irname) == 0 and nelec_fix != mol.nelectron:
                msg =('Num electrons defined by irrep_nelec != total num electrons. '
                      'mol.nelectron = %d  irrep_nelec = %s' %
                      (mol.nelectron, self.irrep_nelec))
                raise ValueError(msg)
            else:
                logger.info(mol, '    %d free electrons in irreps %s',
                            mol.nelectron-nelec_fix, ' '.join(float_irname))
        return ghf.GHF.build(self, mol)

 def shielding(self, mo1=None):
     if getattr(self._scf, 'spin_square', None):
         s2 = self._scf.spin_square()[0]
         if s2 > 1e-4:
             logger.warn(self, '<S^2> = %s. UHF-NMR shielding may have large error.\n'
                         'paramagnetic NMR should include this result plus '
                         'g-tensor and HFC tensors.', s2)
     return rhf_nmr.NMR.shielding(self, mo1)

 def dump_flags(self):
     rhf_mag.Magnetizability.dump_flags(self)
     if self.gauge_orig is not None:
         logger.warn(self, 'Rotational g-tensor with '
                     'perturbation-independent basis is in testing.\n'
                     'Results do not fully agree with those in '
                     'JCP, 105, 2804.')
     return self

 def check_sanity(self):
     mol_hf.SCF.check_sanity(self)
     self.with_df.check_sanity()
     if (isinstance(self.exxdiv, str) and self.exxdiv.lower() != 'ewald' and
         isinstance(self.with_df, df.df.DF)):
         logger.warn(self, 'exxdiv %s is not supported in DF or MDF',
                     self.exxdiv)
     return self

 def dump_flags(self):
     if hasattr(self, "nelectron_alpha"):
         logger.warn(self, "Note the API updates: attribute nelectron_alpha was replaced by attribute nelec")
         # raise RuntimeError('API updates')
         self.nelec = (self.nelectron_alpha, self.mol.nelectron - self.nelectron_alpha)
         delattr(self, "nelectron_alpha")
     hf.SCF.dump_flags(self)
     logger.info(self, "number electrons alpha = %d  beta = %d", *self.nelec)

 def dump_flags(self, verbose=None):
     oldCAS.dump_flags(self)
     self.with_solvent.check_sanity()
     self.with_solvent.dump_flags()
     if self.conv_tol < 1e-7:
         logger.warn(self, 'CASSCF+ddCOSMO may not be able to '
                     'converge to conv_tol=%g', self.conv_tol)
     return self