def shielding(self, mo1=None):
        cput0 = (time.clock(), time.time())
        self.check_sanity()
        self.dump_flags()

        unit_ppm = nist.ALPHA**2 * 1e6
        msc_dia = self.dia(self.gauge_orig)

        if mo1 is None:
            self.mo10, self.mo_e10 = self.solve_mo1()
            mo1 = self.mo10
        msc_para, para_vir, para_occ = self.para(mo10=mo1)

        msc_dia *= unit_ppm
        msc_para *= unit_ppm
        para_vir *= unit_ppm
        para_occ *= unit_ppm
        e11 = msc_para + msc_dia

        logger.timer(self, 'NMR shielding', *cput0)
        if self.verbose >= logger.NOTE:
            for i, atm_id in enumerate(self.shielding_nuc):
                _write(self.stdout, e11[i],
                       '\ntotal shielding of atom %d %s' \
                       % (atm_id, self.mol.atom_symbol(atm_id)))
                _write(self.stdout, msc_dia[i], 'dia-magnetic contribution')
                _write(self.stdout, msc_para[i], 'para-magnetic contribution')
                if self.verbose >= logger.INFO:
                    _write(self.stdout, para_occ[i], 'occ part of para-magnetism')
                    _write(self.stdout, para_vir[i], 'vir part of para-magnetism')
        return e11

    def scf(self, dm0=None):
        cput0 = (time.clock(), time.time())
        mol = self.mol
        self.build(mol)
        self.dump_flags()
        self.converged, self.hf_energy, \
                self.mo_energy, self.mo_coeff, self.mo_occ \
                = hf.kernel(self, self.conv_tol, dm0=dm0,
                            callback=self.callback)

        logger.timer(self, 'SCF', *cput0)
        self.dump_energy(self.hf_energy, self.converged)

        # sort MOs wrt orbital energies, it should be done last.
        o_sort = numpy.argsort(self.mo_energy[self.mo_occ>0])
        v_sort = numpy.argsort(self.mo_energy[self.mo_occ==0])
        self.mo_energy = numpy.hstack((self.mo_energy[self.mo_occ>0][o_sort], \
                                       self.mo_energy[self.mo_occ==0][v_sort]))
        self.mo_coeff = numpy.hstack((self.mo_coeff[:,self.mo_occ>0][:,o_sort], \
                                      self.mo_coeff[:,self.mo_occ==0][:,v_sort]))
        nocc = len(o_sort)
        self.mo_occ[:nocc] = self.mo_occ[self.mo_occ>0][o_sort]
        self.mo_occ[nocc:] = 0

        #if self.verbose >= logger.INFO:
        #    self.analyze(self.verbose)
        return self.hf_energy

    def get_jk(self, mol=None, dm=None, hermi=1):
        if mol is None: mol = self.mol
        if dm is None: dm = self.make_rdm1()
        t0 = (time.clock(), time.time())
        verbose_bak, mol.verbose = mol.verbose, self.verbose
        stdout_bak,  mol.stdout  = mol.stdout , self.stdout
        if self.direct_scf and self.opt[0] is None:
            self.opt = self.init_direct_scf(mol)
        opt_llll, opt_ssll, opt_ssss, opt_gaunt = self.opt

        vj, vk = get_jk_coulomb(mol, dm, hermi, self._coulomb_now,
                                opt_llll, opt_ssll, opt_ssss)

        if self.with_breit:
            if 'SSSS' in self._coulomb_now.upper() or not self.with_ssss:
                vj1, vk1 = _call_veff_gaunt_breit(mol, dm, hermi, opt_gaunt, True)
                logger.info(self, 'Add Breit term')
                vj += vj1
                vk += vk1
        elif self.with_gaunt and 'SS' in self._coulomb_now.upper():
            logger.info(self, 'Add Gaunt term')
            vj1, vk1 = _call_veff_gaunt_breit(mol, dm, hermi, opt_gaunt, False)
            vj += vj1
            vk += vk1

        mol.verbose = verbose_bak
        mol.stdout  = stdout_bak
        logger.timer(self, 'vj and vk', *t0)
        return vj, vk

    def scf(self, dm0=None):
        '''main routine for SCF

        Kwargs:
            dm0 : ndarray
                If given, it will be used as the initial guess density matrix

        Examples:

        >>> import numpy
        >>> from pyscf import gto, scf
        >>> mol = gto.M(atom='H 0 0 0; F 0 0 1.1')
        >>> mf = scf.hf.SCF(mol)
        >>> dm_guess = numpy.eye(mol.nao_nr())
        >>> mf.kernel(dm_guess)
        converged SCF energy = -98.5521904482821
        -98.552190448282104
        '''
        cput0 = (time.clock(), time.time())

        self.build()
        self.dump_flags()
        self.converged, self.hf_energy, \
                self.mo_energy, self.mo_coeff, self.mo_occ = \
                kernel(self, self.conv_tol, dm0=dm0, callback=self.callback)

        logger.timer(self, 'SCF', *cput0)
        self.dump_energy(self.hf_energy, self.converged)
        #if self.verbose >= logger.INFO:
        #    self.analyze(self.verbose)
        return self.hf_energy

def solve_mo1(sscobj, mo_energy=None, mo_coeff=None, mo_occ=None,
              h1=None, s1=None, with_cphf=None):
    cput1 = (time.clock(), time.time())
    log = logger.Logger(sscobj.stdout, sscobj.verbose)
    if mo_energy is None: mo_energy = sscobj._scf.mo_energy
    if mo_coeff  is None: mo_coeff = sscobj._scf.mo_coeff
    if mo_occ    is None: mo_occ = sscobj._scf.mo_occ
    if with_cphf is None: with_cphf = sscobj.cphf

    mol = sscobj.mol
    if h1 is None:
        atmlst = sorted(set([j for i,j in sscobj.nuc_pair]))
        h1 = numpy.asarray(make_h1_pso(mol, mo_coeff, mo_occ, atmlst))

    if with_cphf:
        if callable(with_cphf):
            vind = with_cphf
        else:
            vind = gen_vind(sscobj._scf, mo_coeff, mo_occ)
        mo1, mo_e1 = cphf.solve(vind, mo_energy, mo_occ, h1, None,
                                sscobj.max_cycle_cphf, sscobj.conv_tol,
                                verbose=log)
    else:
        e_ai = lib.direct_sum('i-a->ai', mo_energy[mo_occ>0], mo_energy[mo_occ==0])
        mo1 = h1 * (1 / e_ai)
        mo_e1 = None
    logger.timer(sscobj, 'solving mo1 eqn', *cput1)
    return mo1, mo_e1

    def init_amps(self, eris):
        time0 = time.clock(), time.time()
        nocc = self.nocc()
        nvir = self.nmo() - nocc
        nkpts = self.nkpts
        t1 = numpy.zeros((nkpts,nocc,nvir), dtype=numpy.complex128)
        t2 = numpy.zeros((nkpts,nkpts,nkpts,nocc,nocc,nvir,nvir), dtype=numpy.complex128)
        self.emp2 = 0
        foo = eris.fock[:,:nocc,:nocc].copy()
        fvv = eris.fock[:,nocc:,nocc:].copy()
        eris_oovv = eris.oovv.copy()
        eia = numpy.zeros((nocc,nvir))
        eijab = numpy.zeros((nocc,nocc,nvir,nvir))

        kconserv = tools.get_kconserv(self._scf.cell,self.kpts)
        for ki in range(nkpts):
          for kj in range(nkpts):
            for ka in range(nkpts):
                kb = kconserv[ki,ka,kj]
                for i in range(nocc):
                    for a in range(nvir):
                        eia[i,a] = foo[ki,i,i] - fvv[ka,a,a]
                        for j in range(nocc):
                            for b in range(nvir):
                                eijab[i,j,a,b] = ( foo[ki,i,i] + foo[kj,j,j]
                                                 - fvv[ka,a,a] - fvv[kb,b,b] )
                                t2[ki,kj,ka,i,j,a,b] = eris_oovv[ki,kj,ka,i,j,a,b]/eijab[i,j,a,b]

        t2 = numpy.conj(t2)
        self.emp2 = 0.25*numpy.einsum('pqrijab,pqrijab',t2,eris_oovv).real
        self.emp2 /= nkpts
        logger.info(self, 'Init t2, MP2 energy = %.15g', self.emp2.real)
        logger.timer(self, 'init mp2', *time0)
        print "MP2 energy =", self.emp2
        return self.emp2, t1, t2

    def shielding(self, mo1=None):
        cput0 = (time.clock(), time.time())
        self.dump_flags()
        if self.verbose >= logger.WARN:
            self.check_sanity()

        facppm = 1e6/param.LIGHTSPEED**2
        t0 = (time.clock(), time.time())
        msc_dia = self.dia() * facppm
        t0 = logger.timer(self, 'h11', *t0)
        msc_para, para_pos, para_neg, para_occ = \
                [x*facppm for x in self.para(mo10=mo1)]
        e11 = msc_para + msc_dia

        logger.timer(self, 'NMR shielding', *cput0)
        if self.verbose > param.VERBOSE_QUIET:
            for i, atm_id in enumerate(self.shielding_nuc):
                rhf_nmr._write(self.stdout, e11[i],
                               '\ntotal shielding of atom %d %s'
                               % (atm_id, self.mol.atom_symbol(atm_id-1)))
                rhf_nmr._write(self.stdout, msc_dia[i], 'dia-magnetism')
                rhf_nmr._write(self.stdout, msc_para[i], 'para-magnetism')
                if self.verbose >= param.VERBOSE_INFO:
                    rhf_nmr._write(self.stdout, para_occ[i], 'occ part of para-magnetism')
                    rhf_nmr._write(self.stdout, para_pos[i], 'vir-pos part of para-magnetism')
                    rhf_nmr._write(self.stdout, para_neg[i], 'vir-neg part of para-magnetism')
            self.stdout.flush()
        return e11

    def solve_mo1(self, mo_energy=None, mo_occ=None, h1=None, s1=None):
        cput1 = (time.clock(), time.time())
        log = logger.Logger(self.stdout, self.verbose)
        if mo_energy is None:
            mo_energy = self._scf.mo_energy
        if mo_occ is None:
            mo_occ = self._scf.mo_occ

        mol = self.mol
        if h1 is None:
            mo_coeff = self._scf.mo_coeff
            dm0 = self._scf.make_rdm1(mo_coeff, mo_occ)
            h1 = _mat_ao2mo(self.make_h10(mol, dm0), mo_coeff, mo_occ)
        if s1 is None:
            s1 = _mat_ao2mo(self.make_s10(mol), mo_coeff, mo_occ)

        cput1 = log.timer("first order Fock matrix", *cput1)
        if self.cphf:
            mo10, mo_e10 = cphf.solve(
                self.get_vind, mo_energy, mo_occ, h1, s1, self.max_cycle_cphf, self.conv_tol, verbose=log
            )
        else:
            mo10, mo_e10 = solve_mo1(mo_energy, mo_occ, h1, s1)
        logger.timer(self, "solving mo1 eqn", *cput1)
        return mo10, mo_e10

    def get_jk(self, cell=None, dm=None, hermi=1, kpt=None, kpt_band=None):
        '''Get Coulomb (J) and exchange (K) following :func:`scf.hf.RHF.get_jk_`.

        Note the incore version, which initializes an _eri array in memory.
        '''
        if cell is None: cell = self.cell
        if dm is None: dm = self.make_rdm1()
        if kpt is None: kpt = self.kpt

        cpu0 = (time.clock(), time.time())

        if (kpt_band is None and
            (self.exxdiv == 'ewald' or self.exxdiv is None) and
            (self._eri is not None or cell.incore_anyway or self._is_mem_enough())):
            if self._eri is None:
                logger.debug(self, 'Building PBC AO integrals incore')
                self._eri = self.with_df.get_ao_eri(kpt, compact=True)
            vj, vk = dot_eri_dm(self._eri, dm, hermi)

            if self.exxdiv == 'ewald':
                from pyscf.pbc.df.df_jk import _ewald_exxdiv_for_G0
                # G=0 is not inculded in the ._eri integrals
                _ewald_exxdiv_for_G0(self.cell, kpt, [dm], [vk])
        else:
            vj, vk = self.with_df.get_jk(dm, hermi, kpt, kpt_band,
                                         exxdiv=self.exxdiv)

        logger.timer(self, 'vj and vk', *cpu0)
        return vj, vk

 def get_jk(self, mol=None, dm=None, hermi=1):
     if mol is None: mol = self.mol
     if dm is None: dm = self.make_rdm1()
     t0 = (time.clock(), time.time())
     vj, vk = get_jk(mol, dm, hermi, self.opt)
     logger.timer(self, 'vj and vk', *t0)
     return vj, vk

def grad_elec(mfg, mo_energy=None, mo_coeff=None, mo_occ=None):
    t0 = (time.clock(), time.time())
    mf = mfg._scf
    mol = mfg.mol
    if mo_energy is None: mo_energy = mf.mo_energy
    if mo_occ is None:    mo_occ = mf.mo_occ
    if mo_coeff is None:  mo_coeff = mf.mo_coeff
    h1 = mfg.get_hcore(mol)
    s1 = mfg.get_ovlp(mol)
    dm0 = mf.make_rdm1(mf.mo_coeff, mf.mo_occ)
    vhf = mfg.get_veff(mol, dm0)
    log.timer(mfg, 'gradients of 2e part', *t0)
    f1 = h1 + vhf
    dme0 = mfg.make_rdm1e(mf.mo_energy, mf.mo_coeff, mf.mo_occ)
    gs = numpy.empty((mol.natm,3))
    for ia in range(mol.natm):
# h1, s1, vhf are \nabla <i|h|j>, the nuclear gradients = -\nabla
        f =-(mfg.matblock_by_atom(mol, ia, f1) + mfg._grad_rinv(mol, ia))
        s = -mfg.matblock_by_atom(mol, ia, s1)
        v = numpy.einsum('ij,kji->k', dm0, f) \
          - numpy.einsum('ij,kji->k', dme0, s)
        gs[ia] = 2 * v.real
    log.debug(mfg, 'gradients of electronic part')
    log.debug(mfg, str(gs))
    return gs

    def get_j(self, cell=None, dm=None, hermi=1, kpt=None, kpts_band=None):
        r'''Compute J matrix for the given density matrix and k-point (kpt).
        When kpts_band is given, the J matrices on kpts_band are evaluated.

            J_{pq} = \sum_{rs} (pq|rs) dm[s,r]

        where r,s are orbitals on kpt. p and q are orbitals on kpts_band
        if kpts_band is given otherwise p and q are orbitals on kpt.
        '''
        #return self.get_jk(cell, dm, hermi, kpt, kpts_band)[0]
        if cell is None: cell = self.cell
        if dm is None: dm = self.make_rdm1()
        if kpt is None: kpt = self.kpt

        cpu0 = (time.clock(), time.time())
        dm = np.asarray(dm)
        nao = dm.shape[-1]

        if (kpts_band is None and
            (self._eri is not None or cell.incore_anyway or
             (not self.direct_scf and self._is_mem_enough()))):
            if self._eri is None:
                logger.debug(self, 'Building PBC AO integrals incore')
                self._eri = self.with_df.get_ao_eri(kpt, compact=True)
            vj, vk = mol_hf.dot_eri_dm(self._eri, dm.reshape(-1,nao,nao), hermi)
        else:
            vj = self.with_df.get_jk(dm.reshape(-1,nao,nao), hermi,
                                     kpt, kpts_band, with_k=False)[0]
        logger.timer(self, 'vj', *cpu0)
        return _format_jks(vj, dm, kpts_band)

 def init_amps(self, eris):
     time0 = time.clock(), time.time()
     mo_e = eris.fock.diagonal()
     nocc = self.nocc()
     nvir = mo_e.size - nocc
     t1 = np.zeros((nocc,nvir), eris.dtype)
     #eia = mo_e[:nocc,None] - mo_e[None,nocc:]
     #t1 = eris.fock[:nocc,nocc:] / eia
     t2 = np.zeros((nocc,nocc,nvir,nvir), eris.dtype)
     self.emp2 = 0
     foo = eris.fock[:nocc,:nocc]
     fvv = eris.fock[nocc:,nocc:]
     eia = np.zeros((nocc,nvir))
     eijab = np.zeros((nocc,nocc,nvir,nvir))
     for i in range(nocc):
         for a in range(nvir):
             eia[i,a] = (foo[i,i] - fvv[a,a]).real
             for j in range(nocc):
                 for b in range(nvir):
                     eijab[i,j,a,b] = ( foo[i,i] + foo[j,j]
                                      - fvv[a,a] - fvv[b,b] ).real
                     t2[i,j,a,b] = eris.oovv[i,j,a,b]/eijab[i,j,a,b]
     eris_oovv = _cp(eris.oovv)
     self.emp2 = 0.25*einsum('ijab,ijab',t2,eris_oovv.conj()).real
     logger.info(self, 'Init t2, MP2 energy = %.15g', self.emp2)
     logger.timer(self, 'init mp2', *time0)
     return self.emp2, t1, t2

    def shielding(self, mo1=None):
        cput0 = (time.clock(), time.time())
        self.dump_flags()
        if self.verbose >= logger.WARN:
            self.check_sanity()

        t0 = (time.clock(), time.time())
        unit_ppm = nist.ALPHA**2 * 1e6
        msc_dia = self.dia() * unit_ppm
        t0 = logger.timer(self, 'h11', *t0)
        msc_para, para_pos, para_neg, para_occ = \
                [x*unit_ppm for x in self.para(mo10=mo1)]
        e11 = msc_para + msc_dia

        logger.timer(self, 'NMR shielding', *cput0)
        if self.verbose > logger.QUIET:
            for i, atm_id in enumerate(self.shielding_nuc):
                rhf_nmr._write(self.stdout, e11[i],
                               '\ntotal shielding of atom %d %s'
                               % (atm_id, self.mol.atom_symbol(atm_id)))
                rhf_nmr._write(self.stdout, msc_dia[i], 'dia-magnetism')
                rhf_nmr._write(self.stdout, msc_para[i], 'para-magnetism')
                if self.verbose >= logger.INFO:
                    rhf_nmr._write(self.stdout, para_occ[i], 'occ part of para-magnetism')
                    rhf_nmr._write(self.stdout, para_pos[i], 'vir-pos part of para-magnetism')
                    rhf_nmr._write(self.stdout, para_neg[i], 'vir-neg part of para-magnetism')
        return e11

    def get_jk_(self, cell=None, dm=None, hermi=1, kpt=None, kpt_band=None):
        '''Get Coulomb (J) and exchange (K) following :func:`scf.hf.RHF.get_jk_`.

        Note the incore version, which initializes an _eri array in memory.
        '''
        if cell is None: cell = self.cell
        if dm is None: dm = self.make_rdm1()
        if kpt is None: kpt = self.kpt

        cpu0 = (time.clock(), time.time())
        vj, vk = get_jk(self, cell, dm, hermi, self.opt, kpt, kpt_band)
        # TODO: Check incore, direct_scf, _eri's, etc
        #if self._eri is not None or cell.incore_anyway or self._is_mem_enough():
        #    print "self._is_mem_enough() =", self._is_mem_enough()
        #    if self._eri is None:
        #        logger.debug(self, 'Building PBC AO integrals incore')
        #        if kpt is not None and pyscf.lib.norm(kpt) > 1.e-15:
        #            raise RuntimeError("Non-zero kpts not implemented for incore eris")
        #        self._eri = ao2mo.get_ao_eri(cell)
        #    if np.iscomplexobj(dm) or np.iscomplexobj(self._eri):
        #        vj, vk = dot_eri_dm_complex(self._eri, dm, hermi)
        #    else:
        #        vj, vk = pyscf.scf.hf.dot_eri_dm(self._eri, dm, hermi)
        #else:
        #    if self.direct_scf:
        #        self.opt = self.init_direct_scf(cell)
        #    vj, vk = get_jk(cell, dm, hermi, self.opt, kpt)
        logger.timer(self, 'vj and vk', *cpu0)
        return vj, vk

        def kernel(self, mo_coeff=None, mo_occ=None):
            if mo_coeff is None:
                mo_coeff = self.mo_coeff
            if mo_occ is None:
                mo_occ = self.mo_occ
            cput0 = (time.clock(), time.time())

            self.build(self.mol)
            self.dump_flags()

            if mo_coeff is None or mo_occ is None:
                logger.debug(self, 'Initial guess orbitals not given. '
                             'Generating initial guess from %s density matrix',
                             self.init_guess)
                dm = mf.get_init_guess(self.mol, self.init_guess)
                mo_coeff, mo_occ = self.from_dm(dm)
            # save initial guess because some methods may access them
            self.mo_coeff = mo_coeff
            self.mo_occ = mo_occ

            self.converged, self.e_tot, \
                    self.mo_energy, self.mo_coeff, self.mo_occ = \
                    kernel(self, mo_coeff, mo_occ, conv_tol=self.conv_tol,
                           conv_tol_grad=self.conv_tol_grad,
                           max_cycle=self.max_cycle,
                           callback=self.callback, verbose=self.verbose)

            logger.timer(self, 'Second order SCF', *cput0)
            self._finalize()
            return self.e_tot

    def init_amps(self, eris=None):
        time0 = time.clock(), time.time()
        if eris is None:
            eris = self.ao2mo(self.mo_coeff)
        nocca, noccb = self.nocc

        fova = eris.focka[:nocca,nocca:]
        fovb = eris.fockb[:noccb,noccb:]
        mo_ea_o = eris.mo_energy[0][:nocca]
        mo_ea_v = eris.mo_energy[0][nocca:]
        mo_eb_o = eris.mo_energy[1][:noccb]
        mo_eb_v = eris.mo_energy[1][noccb:]
        eia_a = lib.direct_sum('i-a->ia', mo_ea_o, mo_ea_v)
        eia_b = lib.direct_sum('i-a->ia', mo_eb_o, mo_eb_v)

        t1a = fova.conj() / eia_a
        t1b = fovb.conj() / eia_b

        eris_ovov = np.asarray(eris.ovov)
        eris_OVOV = np.asarray(eris.OVOV)
        eris_ovOV = np.asarray(eris.ovOV)
        t2aa = eris_ovov.transpose(0,2,1,3) / lib.direct_sum('ia+jb->ijab', eia_a, eia_a)
        t2ab = eris_ovOV.transpose(0,2,1,3) / lib.direct_sum('ia+jb->ijab', eia_a, eia_b)
        t2bb = eris_OVOV.transpose(0,2,1,3) / lib.direct_sum('ia+jb->ijab', eia_b, eia_b)
        t2aa = t2aa - t2aa.transpose(0,1,3,2)
        t2bb = t2bb - t2bb.transpose(0,1,3,2)
        e  =      np.einsum('iJaB,iaJB', t2ab, eris_ovOV)
        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)
        self.emp2 = e.real
        logger.info(self, 'Init t2, MP2 energy = %.15g', self.emp2)
        logger.timer(self, 'init mp2', *time0)
        return self.emp2, (t1a,t1b), (t2aa,t2ab,t2bb)

 def get_jk(self, mol=None, dm=None, hermi=0):
     if mol is None: mol = self.mol
     if dm is None: dm = self._scf.make_rdm1()
     cpu0 = (time.clock(), time.time())
     #TODO: direct_scf opt
     vj, vk = get_jk(mol, dm)
     logger.timer(self, 'vj and vk', *cpu0)
     return vj, vk

 def get_jk_(self, mol=None, dm=None, hermi=1):
     if mol is None: mol = self.mol
     if dm is None: dm = self.make_rdm1()
     t0 = (time.clock(), time.time())
     if self.direct_scf and self.opt is None:
         self.opt = self.init_direct_scf(mol)
     vj, vk = get_jk(mol, dm, hermi, self.opt)
     logger.timer(self, 'vj and vk', *t0)
     return vj, vk

def write_chk(mc, root, chkfile):

    t0 = (time.clock(), time.time())
    fh5 = h5py.File(chkfile, "w")

    if mc.fcisolver.nroots > 1:
        mc.mo_coeff, _, mc.mo_energy = mc.canonicalize(mc.mo_coeff, ci=root)

    fh5["mol"] = format(mc.mol.pack())
    fh5["mc/mo"] = mc.mo_coeff
    fh5["mc/ncore"] = mc.ncore
    fh5["mc/ncas"] = mc.ncas
    nvirt = mc.mo_coeff.shape[1] - mc.ncas - mc.ncore
    fh5["mc/nvirt"] = nvirt
    fh5["mc/nelecas"] = mc.nelecas
    fh5["mc/root"] = root
    fh5["mc/orbe"] = mc.mo_energy
    if hasattr(mc, "orbsym"):
        fh5.create_dataset("mc/orbsym", data=mc.orbsym)
    else:
        fh5.create_dataset("mc/orbsym", data=[])

    mo_core = mc.mo_coeff[:, : mc.ncore]
    mo_cas = mc.mo_coeff[:, mc.ncore : mc.ncore + mc.ncas]
    mo_virt = mc.mo_coeff[:, mc.ncore + mc.ncas :]
    core_dm = numpy.dot(mo_core, mo_core.T) * 2
    core_vhf = mc.get_veff(mc.mol, core_dm)
    h1e_Sr = reduce(numpy.dot, (mo_virt.T, mc.get_hcore() + core_vhf, mo_cas))
    h1e_Si = reduce(numpy.dot, (mo_cas.T, mc.get_hcore() + core_vhf, mo_core))
    fh5["h1e_Si"] = h1e_Si
    fh5["h1e_Sr"] = h1e_Sr
    h1e = mc.h1e_for_cas()
    fh5["h1e"] = h1e[0]

    if mc._scf._eri is None:
        from pyscf.scf import _vhf

        eri = _vhf.int2e_sph(mc.mol._atm, mol._bas, mol._env)
    else:
        eri = mc._scf._eri

    # FIXME
    # add outcore later

    h2e = ao2mo.incore.general(eri, [mo_cas, mo_cas, mo_cas, mo_cas], compact=False)
    h2e = h2e.reshape(mc.ncas, mc.ncas, mc.ncas, mc.ncas)
    fh5["h2e"] = h2e
    h2e_Sr = ao2mo.incore.general(eri, [mo_virt, mo_cas, mo_cas, mo_cas], compact=False)
    h2e_Sr = h2e_Sr.reshape(nvirt, mc.ncas, mc.ncas, mc.ncas)
    fh5["h2e_Sr"] = h2e_Sr
    h2e_Si = ao2mo.incore.general(eri, [mo_cas, mo_core, mo_cas, mo_cas], compact=False)
    h2e_Si = h2e_Si.reshape(mc.ncas, mc.ncore, mc.ncas, mc.ncas)
    fh5["h2e_Si"] = h2e_Si

    fh5.close()

    logger.timer(mc, "Write MPS NEVPT integral", *t0)