def make_modchg_basis(auxcell, smooth_eta, l_max=3):
# * chgcell defines smooth gaussian functions for each angular momentum for
#   auxcell. The smooth functions may be used to carry the charge
    chgcell = copy.copy(auxcell)  # smooth model density for coulomb integral to carry charge
    half_sph_norm = .5/numpy.sqrt(numpy.pi)
    chg_bas = []
    chg_env = [smooth_eta]
    ptr_eta = auxcell._env.size
    ptr = ptr_eta + 1
    for ia in range(auxcell.natm):
        for l in set(auxcell._bas[auxcell._bas[:,gto.ATOM_OF]==ia, gto.ANG_OF]):
            if l <= l_max:
                norm = half_sph_norm/gto.mole._gaussian_int(l*2+2, smooth_eta)
                chg_bas.append([ia, l, 1, 1, 0, ptr_eta, ptr, 0])
                chg_env.append(norm)
                ptr += 1

    chgcell._atm = auxcell._atm
    chgcell._bas = numpy.asarray(chg_bas, dtype=numpy.int32).reshape(-1,gto.BAS_SLOTS)
    chgcell._env = numpy.hstack((auxcell._env, chg_env))
    chgcell.nimgs = auxcell.nimgs
    chgcell._built = True
    logger.debug(auxcell, 'smooth basis, num shells = %d, num cGTO = %d',
                 chgcell.nbas, chgcell.nao_nr())
    return chgcell

def kernel(mp, mo_energy, mo_coeff, nocc, ioblk=256, verbose=None):
    nmo = mo_coeff.shape[1]
    nvir = nmo - nocc
    auxmol = df.incore.format_aux_basis(mp.mol, mp.auxbasis)
    naoaux = auxmol.nao_nr()

    iolen = max(int(ioblk*1e6/8/(nvir*nocc)), 160)

    eia = lib.direct_sum('i-a->ia', mo_energy[:nocc], mo_energy[nocc:])
    t2 = None
    emp2 = 0
    with mp.ao2mo(mo_coeff, nocc) as fov:
        for p0, p1 in prange(0, naoaux, iolen):
            logger.debug(mp, 'Load cderi block %d:%d', p0, p1)
            qov = numpy.array(fov[p0:p1], copy=False)
            for i in range(nocc):
                buf = numpy.dot(qov[:,i*nvir:(i+1)*nvir].T,
                                qov).reshape(nvir,nocc,nvir)
                gi = numpy.array(buf, copy=False)
                gi = gi.reshape(nvir,nocc,nvir).transpose(1,2,0)
                t2i = gi/lib.direct_sum('jb+a->jba', eia, eia[i])
                # 2*ijab-ijba
                theta = gi*2 - gi.transpose(0,2,1)
                emp2 += numpy.einsum('jab,jab', t2i, theta)

    return emp2, t2

    def fupdate(t1, t2, istep, normt, de, adiis):
        nocc, nvir = t1.shape
        nov = nocc*nvir
        moidx = numpy.ones(mycc.mo_energy.size, dtype=numpy.bool)
        if isinstance(mycc.frozen, (int, numpy.integer)):
            moidx[:mycc.frozen] = False
        else:
            moidx[mycc.frozen] = False
        mo_e = mycc.mo_energy[moidx]
        eia = mo_e[:nocc,None] - mo_e[None,nocc:]
        if (istep > mycc.diis_start_cycle and
            abs(de) < mycc.diis_start_energy_diff):
            if mycc.t1 is None:
                mycc.t1 = t1
                mycc.t2 = t2
            else:
                tbuf = numpy.empty(nov*(nov+1))
                tbuf[:nov] = ((t1-mycc.t1)*eia).ravel()
                pbuf = tbuf[nov:].reshape(nocc,nocc,nvir,nvir)
                for i in range(nocc):
                    djba = (eia.reshape(-1,1) + eia[i].reshape(1,-1)).reshape(-1)
                    pbuf[i] = (t2[i]-mycc.t2[i]) * djba.reshape(nocc,nvir,nvir)
                adiis.push_err_vec(tbuf)
                tbuf = numpy.empty(nov*(nov+1))
                tbuf[:nov] = t1.ravel()
                tbuf[nov:] = t2.ravel()
                t1.data = tbuf.data # release memory
                t2.data = tbuf.data

                tbuf = adiis.update(tbuf)
                mycc.t1 = t1 = tbuf[:nov].reshape(nocc,nvir)
                mycc.t2 = t2 = tbuf[nov:].reshape(nocc,nocc,nvir,nvir)
            logger.debug(mycc, 'DIIS for step %d', istep)
        return t1, t2

    def assemble_frag_energy(self, mol):
        e_tot = 0
        nelec = 0
        e_corr = 0

        last_frag = -1
        for m, _, _ in self.all_frags:
            if m != last_frag:
                emb = self.embs[m]
                nimp = len(emb.bas_on_frag)
                _, e2frag, dm1 = \
                        self.solver.run(emb, emb._eri, emb.vfit_ci,
                                        with_1pdm=True, with_e2frag=nimp)
                e_frag, nelec_frag = \
                        self.extract_frag_energy(emb, dm1, e2frag)

                log.debug(self, 'fragment %d FCI-in-HF, frag energy = %.12g, E_corr = %.12g, nelec = %.9g', \
                          m, e_frag, e_frag-emb._ehfinhf, nelec_frag)
            e_corr += e_frag-emb._ehfinhf
            e_tot += e_frag
            nelec += nelec_frag
            last_frag = m
        log.info(self, 'sum(e_frag), e_tot = %.9g, nelec_tot = %.9g', \
                  e_tot, nelec)
        return e_tot, e_corr, nelec

def fit_solver(embsys, fock0, nocc, nimp, dm_ref_alpha):
    #fitp = DmFitObj(fock0, nocc, nimp, dm_ref_alpha, v_V, dm_V)
    def _decompress(vfit):
        idx = numpy.tril_indices(nimp)
        v1 = numpy.zeros((nimp, nimp))
        v1[idx] = vfit
        v1[idx[1],idx[0]] = vfit
        return v1

    ec = [0, 0]
    def diff_dm(vfit):
        f = fock0.copy()
        f[:nimp,:nimp] += _decompress(vfit)
        e, c = scipy.linalg.eigh(f)
        dm0 = numpy.dot(c[:nimp,:nocc], c[:nimp,:nocc].T)
        ddm = dm0 - dm_ref_alpha[:nimp,:nimp]
        ec[:] = (e, c)
        return ddm.flatten()

    def jac_ddm(vfit, *args):
        e, c = ec
        x = mat_v_to_mat_dm1(e, c, nocc, nimp, nimp)
        usymm = symm_trans_mat_for_hermit(nimp)
        nn = usymm.shape[0]
        return numpy.dot(x.reshape(-1,nn), usymm)

    x = scipy.optimize.leastsq(diff_dm, numpy.zeros(nimp*(nimp+1)/2),
                               Dfun=jac_ddm, ftol=1e-8)[0]
    log.debug(embsys, 'ddm %s', diff_dm(x))
    return _decompress(x)

    def project(mfmo, init_mo, ncore, s):
        nocc = ncore + casscf.ncas
        mo0core = init_mo[:,:ncore]
        s1 = reduce(numpy.dot, (mfmo.T, s, mo0core))
        idx = numpy.argsort(numpy.einsum('ij,ij->i', s1, s1))
        logger.debug(casscf, 'Core indices %s', str(numpy.sort(idx[-ncore:])))
        # take HF core
        mocore = mfmo[:,numpy.sort(idx[-ncore:])]

        # take projected CAS space
        mocas = init_mo[:,ncore:nocc] \
              - reduce(numpy.dot, (mocore, mocore.T, s, init_mo[:,ncore:nocc]))
        mocc = lo.orth.vec_lowdin(numpy.hstack((mocore, mocas)), s)

        # remove core and active space from rest
        mou = init_mo[:,nocc:] \
            - reduce(numpy.dot, (mocc, mocc.T, s, init_mo[:,nocc:]))
        mo = lo.orth.vec_lowdin(numpy.hstack((mocc, mou)), s)

        if casscf.verbose >= logger.DEBUG:
            s1 = reduce(numpy.dot, (mo[:,ncore:nocc].T, s, mfmo))
            idx = numpy.argwhere(abs(s1) > 0.4)
            for i,j in idx:
                logger.debug(casscf, 'Init guess <mo-CAS|mo-hf>  %d  %d  %12.8f',
                             ncore+i+1, j+1, s1[i,j])
        return mo

    def __debug_hessian_matvec( self ):

        hessian_analytic = np.zeros( [ self.numVars, self.numVars ], dtype=float )
        
        for cnt in range( self.numVars ):
            vector = np.zeros( [ self.numVars ], dtype=float )
            vector[ cnt ] = 1.0
            hessian_analytic[ :, cnt ] = self.__hessian_matvec( vector )

        original_umatrix = np.array( self.u, copy=True )

        stepsize = 1e-8
        self.__set_gradient()
        gradient_ref = np.array( self.gradient, copy=True )
        hessian_numerical = np.zeros( [ self.numVars, self.numVars ], dtype=float )
        for counter in range( self.numVars ):
            self.u = np.array( original_umatrix, copy=True )
            flatx = np.zeros( [ self.numVars ], dtype=float )
            flatx[counter] = stepsize
            self.__update_unitary( flatx )
            self.__set_gradient()
            hessian_numerical[ :, counter ] = ( self.gradient - gradient_ref ) / stepsize

        self.u = np.array( original_umatrix, copy=True )
        flatx = np.zeros( [ self.numVars ], dtype=float )
        self.__update_unitary( flatx )

        hessian_numerical = 0.5 * ( hessian_numerical + hessian_numerical.T )

        logger.debug(self, "2-norm( hessian difference ) = %g", np.linalg.norm( hessian_analytic - hessian_numerical ))
        logger.debug(self, "2-norm( hessian )            = %g", np.linalg.norm( hessian_analytic ))

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

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

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

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

        def 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 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 get_occ(self, mo_energy=None, mo_coeff=None):
     if mo_energy is None:
         mo_energy = self.mo_energy
     mol = self.mol
     n4c = len(mo_energy)
     n2c = n4c // 2
     mo_occ = numpy.zeros(n2c * 2)
     if mo_energy[n2c] > -1.999 * mol.light_speed ** 2:
         mo_occ[n2c : n2c + mol.nelectron] = 1
     else:
         n = 0
         for i, e in enumerate(mo_energy):
             if e > -1.999 * mol.light_speed ** 2 and n < mol.nelectron:
                 mo_occ[i] = 1
                 n += 1
     if self.verbose >= logger.INFO:
         logger.info(
             self,
             "HOMO %d = %.12g, LUMO %d = %.12g,",
             (n2c + mol.nelectron) // 2,
             mo_energy[n2c + mol.nelectron - 1],
             (n2c + mol.nelectron) // 2 + 1,
             mo_energy[n2c + mol.nelectron],
         )
         logger.debug(self, "NES  mo_energy = %s", mo_energy[:n2c])
         logger.debug(self, "PES  mo_energy = %s", mo_energy[n2c:])
     return mo_occ

def energy_elec(ks, dm, h1e=None, vhf=None):
    if h1e is None:
        h1e = ks.get_hcore()
    e1 = numpy.einsum('ij,ij', h1e.conj(), dm[0]+dm[1])
    tot_e = e1 + ks._ecoul + ks._exc
    logger.debug(ks, 'Ecoul = %s  Exc = %s', ks._ecoul, ks._exc)
    return tot_e, ks._ecoul+ks._exc

def remove_linear_dep_(mf, threshold=LINEAR_DEP_THRESHOLD,
                       lindep=LINEAR_DEP_TRIGGER):
    '''
    Args:
        threshold : float
            The threshold under which the eigenvalues of the overlap matrix are
            discarded to avoid numerical instability.
        lindep : float
            The threshold that triggers the special treatment of the linear
            dependence issue.
    '''
    s = mf.get_ovlp()
    cond = numpy.max(lib.cond(s))
    if cond < 1./lindep:
        return mf

    logger.info(mf, 'Applying remove_linear_dep_ on SCF obejct.')
    logger.debug(mf, 'Overlap condition number %g', cond)
    def eigh(h, s):
        d, t = numpy.linalg.eigh(s)
        x = t[:,d>threshold] / numpy.sqrt(d[d>threshold])
        xhx = reduce(numpy.dot, (x.T.conj(), h, x))
        e, c = numpy.linalg.eigh(xhx)
        c = numpy.dot(x, c)
        return e, c
    mf._eigh = eigh
    return mf

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_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 build_(self):
        mol = self.mol
        self.orth_coeff = self.get_orth_ao(mol)
        self.bas_on_frag = select_ao_on_fragment(mol, self.imp_atoms, \
                                                 self.imp_basidx)
        c_inv = numpy.dot(self.orth_coeff.T, self.entire_scf.get_ovlp(self.mol))
        mo_a, mo_b = self.entire_scf.mo_coeff
        occ_a, occ_b = self.entire_scf.mo_occ
        mo_orth_a = numpy.dot(c_inv, mo_a[:,self.entire_scf.mo_occ[0]>1e-15])
        mo_orth_b = numpy.dot(c_inv, mo_b[:,self.entire_scf.mo_occ[1]>1e-15])
        # self.imp_site, self.bath_orb, self.env_orb are based on orth-orbitals
        self.imp_site, self.bath_orb, self.env_orb = \
                self.decompose_orbital((mo_orth_a,mo_orth_b))
        ovlp = numpy.dot(self.bath_orb[0].T,self.bath_orb[1])[:4,:4]
        for i,c in enumerate(ovlp):
            log.debug(self, ('<bath_alpha_%d|bath_beta> = ' % i) \
                      + '%10.5f'*len(c), *c)
        self.impbas_coeff = self.cons_impurity_basis()
        self.nelectron_alpha = self.entire_scf.nelectron_alpha \
                - self.env_orb[0].shape[1]
        self.nelectron_beta = mol.nelectron \
                - self.entire_scf.nelectron_alpha \
                - self.env_orb[1].shape[1]
        log.info(self, 'alpha / beta electrons = %d / %d', \
                 self.nelectron_alpha, self.nelectron_beta)

        self.energy_by_env, self._vhf_env = self.init_vhf_env(self.env_orb)

    def kernel(self, h1e, eri, norb, nelec, ci0=None, **kwargs):
        if self.verbose > logger.QUIET:
            pyscf.gto.mole.check_sanity(self, self._keys, self.stdout)

        wfnsym = _id_wfnsym(self, norb, nelec, self.wfnsym)
        if 'verbose' in kwargs:
            if isinstance(kwargs['verbose'], logger.Logger):
                log = kwargs['verbose']
            else:
                log = logger.Logger(self.stdout, kwargs['verbose'])
            log.debug('total symmetry = %s',
                      symm.irrep_id2name(self.mol.groupname, wfnsym))
        else:
            logger.debug(self, 'total symmetry = %s',
                         symm.irrep_id2name(self.mol.groupname, wfnsym))
        e, c = direct_spin1.kernel_ms1(self, h1e, eri, norb, nelec, ci0,
                                       **kwargs)
        if self.wfnsym is not None:
            # should I remove the non-symmetric contributions in each
            # call of contract_2e?
            if self.nroots > 1:
                c = [addons.symmetrize_wfn(ci, norb, nelec, self.orbsym, wfnsym)
                     for ci in c]
            else:
                c = addons.symmetrize_wfn(c, norb, nelec, self.orbsym, wfnsym)
        return e, c

 def max_stepsize_scheduler(self, envs):
     if envs['de'] < self.conv_tol or self._max_stepsize is None:
         self._max_stepsize = self.max_stepsize
     else:
         self._max_stepsize *= .5
         logger.debug(self, 'set max_stepsize to %g', self._max_stepsize)
     return self._max_stepsize

    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 vfit_mf_method(self, mol, embsys):
        dm_ref = []
        for m,emb in enumerate(self.embs):
            nimp = len(emb.bas_on_frag)
            dv = numpy.eye(nimp) * embsys.vfit_ci[0]
            dmci = embsys.solver.run(emb, emb._eri, dv, True, False)[2]
            log.debug(embsys, 'dm_ref %d = %s', m, dmci)
            dm_ref.append(dmci*.5)
        if self.translational:
            dm_ref = [dm_ref[0] for i in self.basidx_group]

        sc = reduce(numpy.dot, (self.orth_coeff.T,
                                self.entire_scf.get_ovlp(),
                                self.entire_scf.mo_coeff))
        # this fock matrix includes the previous fitting potential
        fock0 = numpy.dot(sc*self.entire_scf.mo_energy, sc.T.conj())
        nocc = mol.nelectron // 2

        dv = embsys.fit_solver(fock0, nocc, dm_ref)
        v = [embsys.vfit_mf[m]+dv[m] for m,emb in enumerate(self.embs)]
        log.debug(self, 'vfit_mf = ')
        try:
            v_group = self.pack(v)
            for vi in v_group:
                pyscf.tools.dump_mat.dump_tri(self.stdout, vi)
        except:
            self.stdout.write('%s\n' % str(v))
        return v