def _symmetrize_canonicalization_(mol, mo_energy, mo_coeff, s):
    '''Restore symmetry for canonicalized orbitals
    '''
    def search_for_degeneracy(mo_energy):
        idx = numpy.where(abs(mo_energy[1:] - mo_energy[:-1]) < 1e-6)[0]
        return numpy.unique(numpy.hstack((idx, idx+1)))

    degidx = search_for_degeneracy(mo_energy)
    logger.debug1(mol, 'degidx %s', degidx)
    if degidx.size > 0:
        esub = mo_energy[degidx]
        csub = mo_coeff[:,degidx]
        scsub = numpy.dot(s, csub)
        emin = abs(esub).min() * .5
        es = []
        cs = []
        for i,ir in enumerate(mol.irrep_id):
            so = mol.symm_orb[i]
            sosc = numpy.dot(so.T, scsub)
            s_ir = reduce(numpy.dot, (so.T, s, so))
            fock_ir = numpy.dot(sosc*esub, sosc.T)
            mo_energy, u = scipy.linalg.eigh(fock_ir, s_ir)
            idx = abs(mo_energy) > emin
            es.append(mo_energy[idx])
            cs.append(numpy.dot(mol.symm_orb[i], u[:,idx]))
        es = numpy.hstack(es).round(9)
        idx = numpy.argsort(es)
        assert(numpy.allclose(es[idx], esub))
        mo_coeff[:,degidx] = numpy.hstack(cs)[:,idx]
    return mo_coeff

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 dump_flags(self):
     log = logger.Logger(self.stdout, self.verbose)
     logger.info(self, '\n')
     logger.info(self, '******** %s flags ********', self.__class__)
     logger.info(self, 'gs = %s', self.gs)
     logger.info(self, 'len(kpts) = %d', len(self.kpts))
     logger.debug1(self, '    kpts = %s', self.kpts)

 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 aux_e2(cell, auxcell, intor):
    '''3-center AO integrals (ij|L), where L is the auxiliary basis.

    Implements double summation over lattice vectors: \sum_{lm} (i[l]j[m]|L[0]).
    '''
    # sum over largest number of images in either cell or auxcell
    nimgs = numpy.max((cell.nimgs, auxcell.nimgs), axis=0)
    Ls = tools.pbc.get_lattice_Ls(cell, nimgs)
    logger.debug1(cell, "Images summed over in DFT %s", nimgs)
    logger.debug2(cell, "Ls = %s", Ls)

    nao = cell.nao_nr()
    nao_pair = nao*(nao+1) // 2
    nao_pair = nao*nao
    naoaux = auxcell.nao_nr()
    cellL = cell.copy()
    cellR = cell.copy()
    _envL = cellL._env
    _envR = cellR._env
    ptr_coord = cellL._atm[:,pyscf.gto.PTR_COORD]
    buf = numpy.zeros((nao_pair,naoaux))
    for l, L1 in enumerate(Ls):
        _envL[ptr_coord+0] = cell._env[ptr_coord+0] + L1[0]
        _envL[ptr_coord+1] = cell._env[ptr_coord+1] + L1[1]
        _envL[ptr_coord+2] = cell._env[ptr_coord+2] + L1[2]
        for m in range(l):
            _envR[ptr_coord+0] = cell._env[ptr_coord+0] + Ls[m][0]
            _envR[ptr_coord+1] = cell._env[ptr_coord+1] + Ls[m][1]
            _envR[ptr_coord+2] = cell._env[ptr_coord+2] + Ls[m][2]

            buf += pyscf.df.incore.aux_e2(cellL, auxcell, intor, mol1=cellR)
        buf += .5 * pyscf.df.incore.aux_e2(cellL, auxcell, intor, mol1=cellL)
    eri = buf.reshape(nao,nao,-1)
    return eri + eri.transpose(1,0,2)

def make_modchg_basis(auxcell, smooth_eta):
# * 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
    l_max = auxcell._bas[:,gto.ANG_OF].max()
# gaussian_int(l*2+2) for multipole integral:
# \int (r^l e^{-ar^2} * Y_{lm}) (r^l Y_{lm}) r^2 dr d\Omega
    norms = [half_sph_norm/gto.gaussian_int(l*2+2, smooth_eta)
             for l in range(l_max+1)]
    for ia in range(auxcell.natm):
        for l in set(auxcell._bas[auxcell._bas[:,gto.ATOM_OF]==ia, gto.ANG_OF]):
            chg_bas.append([ia, l, 1, 1, 0, ptr_eta, ptr, 0])
            chg_env.append(norms[l])
            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.rcut = _estimate_rcut(smooth_eta, l_max, 1., auxcell.precision)
    logger.debug1(auxcell, 'make compensating basis, num shells = %d, num cGTOs = %d',
                  chgcell.nbas, chgcell.nao_nr())
    logger.debug1(auxcell, 'chgcell.rcut %s', chgcell.rcut)
    return chgcell

 def update(self, s, d, f, *args, **kwargs):
     errvec = get_err_vec(s, d, f)
     logger.debug1(self, 'diis-norm(errvec)=%g', numpy.linalg.norm(errvec))
     xnew = lib.diis.DIIS.update(self, f, xerr=errvec)
     if self.rollback > 0 and len(self._bookkeep) == self.space:
         self._bookkeep = self._bookkeep[-self.rollback:]
     return xnew

def init_guess_by_atom(mol):
    '''Generate initial guess density matrix from superposition of atomic HF
    density matrix.  The atomic HF is occupancy averaged RHF

    Returns:
        Density matrix, 2D ndarray
    '''
    from pyscf.scf import atom_hf
    atm_scf = atom_hf.get_atm_nrhf(mol)
    nbf = mol.nao_nr()
    dm = numpy.zeros((nbf, nbf))
    p0 = 0
    for ia in range(mol.natm):
        symb = mol.atom_symbol(ia)
        if symb in atm_scf:
            e_hf, mo_e, mo_c, mo_occ = atm_scf[symb]
        else:
            symb = mol.atom_pure_symbol(ia)
            e_hf, mo_e, mo_c, mo_occ = atm_scf[symb]
        p1 = p0 + mo_e.__len__()
        dm[p0:p1,p0:p1] = numpy.dot(mo_c*mo_occ, mo_c.T.conj())
        p0 = p1

    for k, v in atm_scf.items():
        logger.debug1(mol, 'Atom %s, E = %.12g', k, v[0])
    return dm

 def dump_flags(self):
     grad_class.dump_flags(self)
     logger.info(self, '** Add background charges for %s **', grad_class)
     if self.verbose >= logger.DEBUG1:
         logger.debug1(self, 'Charge      Location')
         for i, z in enumerate(charges):
             logger.debug1(self, '%.9g    %s', z, coords[i])
     return self

 def dump_flags(self):
     logger.info(self, '\n')
     logger.info(self, '******** %s ********', self.__class__)
     logger.info(self, 'mesh = %s (%d PWs)', self.mesh, numpy.prod(self.mesh))
     logger.info(self, 'eta = %s', self.eta)
     logger.info(self, 'len(kpts) = %d', len(self.kpts))
     logger.debug1(self, '    kpts = %s', self.kpts)
     return self

 def make_rdm1(self, mo_coeff=None, mo_occ=None):
     if mo_coeff is None:
         mo_coeff = self.mo_coeff_on_imp
     if mo_occ is None:
         mo_occ = self.mo_occ
     nbf = mo_coeff.shape[0]
     mo = mo_coeff[:,mo_occ>0]
     dm = numpy.dot(mo, mo.T.conj()) * 2
     log.debug1(self, 'density.diag = %s', dm.diagonal())
     return dm

    def update(self, x, xerr=None):
        '''Extrapolate vector 

        * If xerr the error vector is given, this function will push the target
        vector and error vector in the DIIS subspace, and use the error vector
        to extrapolate the vector and return the extrapolated vector.
        * If xerr is None, this function will take the difference between
        the current given vector and the last given vector as the error
        vector to extrapolate the vector.
        '''
        if xerr is not None:
            self.push_err_vec(xerr)
        self.push_vec(x)

        nd = self.get_num_vec()
        if nd < self.min_space:
            return x

        dt = numpy.array(self.get_err_vec(self._head-1), copy=False)
        for i in range(nd):
            tmp = 0
            dti = self.get_err_vec(i)
            for p0,p1 in prange(0, dt.size, BLOCK_SIZE):
                tmp += numpy.dot(dt[p0:p1].conj(), dti[p0:p1])
            self._H[self._head,i+1] = tmp
            self._H[i+1,self._head] = tmp.conjugate()
        dt = None
        h = self._H[:nd+1,:nd+1]
        g = numpy.zeros(nd+1, x.dtype)
        g[0] = 1

        #try:
        #    c = numpy.linalg.solve(h, g)
        #except numpy.linalg.linalg.LinAlgError:
        #    logger.warn(self, ' diis singular')
        if 1:
            w, v = scipy.linalg.eigh(h)
            idx = abs(w)>1e-14
            c = numpy.dot(v[:,idx]*(1/w[idx]), numpy.dot(v[:,idx].T.conj(), g))
        logger.debug1(self, 'diis-c %s', c)

        if self._xprev is None:
            xnew = numpy.zeros_like(x.ravel())
        else:
            self._xprev = None # release memory first
            self._xprev = xnew = numpy.zeros_like(x.ravel())

        for i, ci in enumerate(c[1:]):
            xi = self.get_vec(i)
            for p0,p1 in prange(0, x.size, BLOCK_SIZE):
                xnew[p0:p1] += xi[p0:p1] * ci
        return xnew.reshape(x.shape)

 def dump_flags(self):
     log = logger.Logger(self.stdout, self.verbose)
     logger.info(self, '\n')
     logger.info(self, '******** %s flags ********', self.__class__)
     logger.info(self, 'gs = %s', self.gs)
     logger.info(self, 'auxbasis = %s', self.auxbasis)
     logger.info(self, 'eta = %s', self.eta)
     if isinstance(self._cderi, str):
         logger.info(self, '_cderi = %s', self._cderi)
     else:
         logger.info(self, '_cderi = %s', self._cderi_file.name)
     logger.info(self, 'len(kpts) = %d', len(self.kpts))
     logger.debug1(self, '    kpts = %s', self.kpts)

def DMRG_MPS_NEVPT(mc, root=0, fcisolver=None,maxm = 500, tol =1e-6, parallel= True):
    
    if (isinstance(mc, basestring)):
        fh5 = h5py.File(mc,'r')

        mol = eval(fh5['mol'].value)
        ncas = fh5['mc/ncas'].value
        ncore = fh5['mc/ncore'].value
        nvirt = fh5['mc/nvirt'].value
        nelecas = fh5['mc/nelecas'].value
        fh5.close()
        mc_chk = mc
    else :
        mol = mc.mol
        ncas = mc.ncas
        ncore = mc.ncore
        nvirt = mc.mo_coeff.shape[1] - mc.ncas-mc.ncore
        nelecas = mc.nelecas
        mc_chk = 'mc_chkfile'
        write_chk(mc,root,mc_chk)
        



    if fcisolver is None:
        fcisolver = DMRGCI(mol, maxm, tol)
    fcisolver.twopdm = False
    fcisolver.nroots = mc.fcisolver.nroots
    scratch = fcisolver.scratchDirectory
    fcisolver.scratchDirectory = ''
    #if (not parallel):
    #    ci.extraline.append('restart_mps_nevpt %d %d %d'%(ncas,ncore, nvirt))


    fcisolver.extraline.append('fullrestart')
    fcisolver.extraline.append('nevpt_state_num %d'%root)
    
    writeDMRGConfFile(nelecas[0], nelecas[1], False, fcisolver)
    fcisolver.scratchDirectory = scratch

    if fcisolver.verbose >= logger.DEBUG1:
        inFile = fcisolver.configFile
        #inFile = os.path.join(self.scratchDirectory,self.configFile)
        logger.debug1(fcisolver, 'Block Input conf')
        logger.debug1(fcisolver, open(inFile, 'r').read())


    from subprocess import check_call
    import os
    full_path = os.path.realpath(__file__)
    check_call('%s %s/nevpt_mpi.py %s %s %s %s %s'%(fcisolver.mpiprefix, os.path.dirname(full_path), mc_chk, fcisolver.executable, fcisolver.configFile,fcisolver.outputFile, fcisolver.scratchDirectory), shell=True)

 def update(self, s, d, f):
     if isinstance(f, numpy.ndarray) and f.ndim == 2:
         sdf = reduce(numpy.dot, (s,d,f))
         errvec = sdf.T.conj() - sdf
     else:
         sdf_a = reduce(numpy.dot, (s, d[0], f[0]))
         sdf_b = reduce(numpy.dot, (s, d[1], f[1]))
         errvec = numpy.hstack((sdf_a.T.conj() - sdf_a,
                                sdf_b.T.conj() - sdf_b))
     logger.debug1(self, 'diis-norm(errvec)=%g', numpy.linalg.norm(errvec))
     xnew = pyscf.lib.diis.DIIS.update(self, f, xerr=errvec)
     if self.rollback > 0 and len(self._bookkeep) == self.space:
         self._bookkeep = self._bookkeep[-self.rollback:]
     return xnew

    def kernel(self, h1e, eri, norb, nelec, ci0=None, ecore=0, restart=None,
               **kwargs):
        if restart is None:
            restart = self.restart
        state_id = min(self.config['eps_vars'])

        if restart or ci0 is not None:
            if self.verbose >= logger.DEBUG1:
                logger.debug1(self, 'restart was set. wf is read from wf_eps* file.')
            self.cleanup(remove_wf=False)
            wfn_file = get_wfn_file(self, state_id)
            if os.path.isfile(wfn_file):
                shutil.move(wfn_file, get_wfn_file(self, state_id * 2))
        else:
            self.cleanup(remove_wf=True)

        if 'orbsym' in kwargs:
            self.orbsym = kwargs['orbsym']
        writeIntegralFile(self, h1e, eri, norb, nelec, ecore)
        conf = {}
        if 'tol' in kwargs:
            conf['tol'] = kwargs['tol']
        write_config(self, nelec, conf)

        if self.dryrun:
            logger.info(self, 'Only write integrals and config')
            if self.nroots == 1:
                calc_e = 0.0
                roots = ''
            else :
                calc_e = [0.0] * self.nroots
                roots = [''] * self.nroots
            return calc_e, roots

        if self.nroots != 1:
            raise NotImplementedError

        execute_shci(self)
        if self.verbose >= logger.DEBUG1:
            with open(os.path.join(self.runtimedir, self.outputfile), 'r') as f:
                self.stdout.write(f.read())

        calc_e = read_energy(self)

        # Each eps_vars is associated to one approximate wfn.
        roots = state_id = min(self.config['eps_vars'])
        if not os.path.isfile(get_wfn_file(self, state_id)):
            raise RuntimeError('Eigenstate %s not found' % get_wfn_file(self, state_id))
        return calc_e, roots

 def dump_flags(self):
     log = logger.Logger(self.stdout, self.verbose)
     logger.info(self, "\n")
     logger.info(self, "******** %s flags ********", self.__class__)
     logger.info(self, "gs = %s", self.gs)
     logger.info(self, "metric = %s", self.metric)
     logger.info(self, "approx_sr_level = %s", self.approx_sr_level)
     logger.info(self, "auxbasis = %s", self.auxbasis)
     logger.info(self, "eta = %s", self.eta)
     if isinstance(self._cderi, str):
         logger.info(self, "_cderi = %s", self._cderi)
     else:
         logger.info(self, "_cderi = %s", self._cderi_file.name)
     logger.info(self, "len(kpts) = %d", len(self.kpts))
     logger.debug1(self, "    kpts = %s", self.kpts)

def get_veff(ks_grad, mol=None, dm=None):
    '''Coulomb + XC functional
    '''
    if mol is None: mol = ks_grad.mol
    if dm is None: dm = ks_grad.base.make_rdm1()
    t0 = (time.clock(), time.time())

    mf = ks_grad.base
    ni = mf._numint
    if ks_grad.grids is not None:
        grids = ks_grad.grids
    else:
        grids = mf.grids
    if grids.coords is None:
        grids.build(with_non0tab=True)

    if mf.nlc != '':
        raise NotImplementedError
    #enabling range-separated hybrids
    omega, alpha, hyb = ni.rsh_and_hybrid_coeff(mf.xc, spin=mol.spin)

    mem_now = lib.current_memory()[0]
    max_memory = max(2000, ks_grad.max_memory*.9-mem_now)
    if ks_grad.grid_response:
        exc, vxc = get_vxc_full_response(ni, mol, grids, mf.xc, dm,
                                         max_memory=max_memory,
                                         verbose=ks_grad.verbose)
        logger.debug1(ks_grad, 'sum(grids response) %s', exc.sum(axis=0))
    else:
        exc, vxc = get_vxc(ni, mol, grids, mf.xc, dm,
                           max_memory=max_memory, verbose=ks_grad.verbose)
    t0 = logger.timer(ks_grad, 'vxc', *t0)

    if abs(hyb) < 1e-10:
        vj = ks_grad.get_j(mol, dm)
        vxc += vj[0] + vj[1]
    else:
        vj, vk = ks_grad.get_jk(mol, dm)
        vk *= hyb
        if abs(omega) > 1e-10:  # For range separated Coulomb operator
            with mol.with_range_coulomb(omega):
                vk += ks_grad.get_k(mol, dm) * (alpha - hyb)
        vxc += vj[0] + vj[1] - vk

    return lib.tag_array(vxc, exc1_grid=exc)

    def eig(self, f, s):
        atm = self.mol
        symb = atm.atom_symbol(0)
        idx_by_l = [[] for i in range(param.L_MAX)]
        i0 = 0
        for ib in range(atm.nbas):
            l = atm.bas_angular(ib)
            nc = atm.bas_nctr(ib)
            i1 = i0 + nc * (l*2+1)
            idx_by_l[l].extend(range(i0, i1, l*2+1))
            i0 = i1

        nbf = atm.nao_nr()
        self._occ = numpy.zeros(nbf)
        mo_c = numpy.zeros((nbf, nbf))
        mo_e = numpy.zeros(nbf)

        # fraction occupation
        for l in range(param.L_MAX):
            if idx_by_l[l]:
                n2occ, frac = frac_occ(symb, l)
                logger.debug1(self, 'l = %d  occ = %d + %.4g', l, n2occ, frac)

                idx = numpy.array(idx_by_l[l])
                f1 = 0
                s1 = 0
                for m in range(l*2+1):
                    f1 = f1 + f[idx+m,:][:,idx+m]
                    s1 = s1 + s[idx+m,:][:,idx+m]
                f1 *= 1./(l*2+1)
                s1 *= 1./(l*2+1)
                e, c = hf.SCF.eig(self, f1, s1)
                for i, ei in enumerate(e):
                    logger.debug1(self, 'l = %d  e_%d = %.9g', l, i, ei)

                for m in range(l*2+1):
                    mo_e[idx] = e
                    self._occ[idx[:n2occ]] = 2
                    if frac > 1e-15:
                        self._occ[idx[n2occ]] = frac
                    for i,i1 in enumerate(idx):
                        mo_c[idx,i1] = c[:,i]
                    idx += 1
        return mo_e, mo_c

    def update(self, x, xerr=None):
        '''use DIIS method to solve Eq.  operator(x) = x.'''
        if xerr is not None:
            self.push_err_vec(xerr)
        self.push_vec(x)

        nd = self.get_num_vec()
        if nd < self.min_space:
            return x

        dt = numpy.array(self.get_err_vec(self._head-1), copy=False)
        for i in range(nd):
            tmp = 0
            dti = self.get_err_vec(i)
            for p0,p1 in prange(0, x.size, BLOCK_SIZE):
                tmp += numpy.dot(dt[p0:p1].conj(), dti[p0:p1])
            self._H[self._head,i+1] = tmp
            self._H[i+1,self._head] = tmp.conjugate()
        dt = None
        h = self._H[:nd+1,:nd+1]
        g = numpy.zeros(nd+1, x.dtype)
        g[0] = 1

        try:
            c = numpy.linalg.solve(h, g)
        except numpy.linalg.linalg.LinAlgError:
            logger.warn(self, 'singularity in diis')
            w, v = scipy.linalg.eigh(h)
            idx = abs(w)>1e-14
            c = numpy.dot(v[:,idx]*(1/w[idx]), numpy.dot(v[:,idx].T.conj(), g))
        logger.debug1(self, 'diis-c %s', c)

        if self._xprev is None:
            xnew = numpy.zeros_like(x.ravel())
        else:
            self._xprev = None # release memory first
            self._xprev = xnew = numpy.zeros_like(x.ravel())

        for i, ci in enumerate(c[1:]):
            xi = self.get_vec(i)
            for p0,p1 in prange(0, x.size, BLOCK_SIZE):
                xnew[p0:p1] += xi[p0:p1] * ci
        return xnew.reshape(x.shape)