def temporal_FT(Qktau, beta, ntau, n_iw, nk, statistic='Fermion', use_IBZ_symmetry = True, N_cores=1):        
  if N_cores==1: return temporal_FT_single_core(Qktau, beta, ntau, n_iw, nk, statistic, use_IBZ_symmetry)
  if mpi.is_master_node(): print "temporal_FT, N_cores: ",N_cores
  if statistic=='Fermion':
    nw = 2*n_iw
  elif statistic=='Boson':
    nw = 2*n_iw-1
  else:
    if mpi.is_master_node(): 
      print "statistic not implemented"
    quit()

  Qkw = numpy.zeros((nw,nk,nk), dtype=numpy.complex_)
        
  if use_IBZ_symmetry: max_kxi = nk/2+1
  else: max_kxi = nk
  pool = Pool(processes=N_cores)              # start worker processes             
  for kxi in range(max_kxi):                    
    if use_IBZ_symmetry: max_kyi = nk/2+1
    else: max_kyi = nk
    numpy.transpose(Qkw[:,kxi,:])[0:max_kyi,:] = pool.map(f_,
                                                      [( deepcopy(Qktau[:,kxi,kyi]),
                                                         beta, ntau, n_iw, statistic
                                                       )\
                                                       for kyi in range(max_kyi)])  
  pool.close()      

  if use_IBZ_symmetry: 
    for wi in range(nw):
      IBZ.copy_by_weak_symmetry(Qkw[wi,:,:], nk)
  return Qkw

    def check_and_fix(self, data, finalize = True, keep_P_negative = True):
      #safe_values = self.get_safe_values(data.Jq, data.bosonic_struct, data.n_q, data.n_q)
      safe_values = {}
      for A in data.bosonic_struct.keys():
        safe_values[A] = 1.0/numpy.amin(data.Jq[A])     
      if mpi.is_master_node(): print "edmft.cautionary: safe_values: ", safe_values
      #print "[Node",mpi.rank,"]","edmft.cautionary: actual safe values: (0,1) = ", 1.0/numpy.amin(data.Jq['0']),1.0/numpy.amin(data.Jq['1'])
      #operates directly on data.P_loc_iw as this is the one that will be used in chiqnu calculation
      clipped = False
  
      prefactor = 1.0 - self.ms0 / (self.clip_counter**self.ccpower + 1.0)
      for A in data.bosonic_struct.keys():
        for i in range(data.nnu):
          if keep_P_negative:
            if (data.P_loc_iw[A].data[i,0,0].real > 0):      
              data.P_loc_iw[A].data[i,0,0] = 0.0
            #clipped = True        
          if (data.P_loc_iw[A].data[i,0,0].real < safe_values[A]) and (safe_values[A]<0.0):      
            data.P_loc_iw[A].data[i,0,0] = prefactor*safe_values[A] + 1j*data.P_loc_iw[A].data[i,0,0].imag
            clipped = True        
            if mpi.is_master_node(): print "edmft.cautionary: clipping P_loc in block ",A
      if clipped and finalize: 
        self.clip_counter += 1 
      else: 
        self.clip_counter = self.clip_counter/self.ccrelax 

      return clipped

    def convert_gf(self,G,G_struct,ish=0,show_warnings=True,**kwargs):
        """ Convert BlockGf from its structure to this structure.

        .. warning::

            Elements that are zero in the new structure due to
            the new block structure will be just ignored, thus
            approximated to zero.

        Parameters
        ----------
        G : BlockGf
            the Gf that should be converted
        G_struct : GfStructure
            the structure of that G
        ish : int
            shell index
        show_warnings : bool or float
            whether to show warnings when elements of the Green's
            function get thrown away
            if float, set the threshold for the magnitude of an element
            about to be thrown away to trigger a warning
            (default: 1.e-10)
        **kwargs :
            options passed to the constructor for the new Gf
        """

        warning_threshold = 1.e-10
        if isinstance(show_warnings, float):
            warning_threshold = show_warnings
            show_warnings = True

        G_new = self.create_gf(ish=ish,**kwargs)
        for block in G_struct.gf_struct_solver[ish].keys():
            for i1 in G_struct.gf_struct_solver[ish][block]:
                for i2 in G_struct.gf_struct_solver[ish][block]:
                    i1_sumk = G_struct.solver_to_sumk[ish][(block,i1)]
                    i2_sumk = G_struct.solver_to_sumk[ish][(block,i2)]
                    i1_sol = self.sumk_to_solver[ish][i1_sumk]
                    i2_sol = self.sumk_to_solver[ish][i2_sumk]
                    if i1_sol[0] is None or i2_sol[0] is None:
                        if show_warnings:
                            if mpi.is_master_node():
                                warn(('Element {},{} of block {} of G is not present '+
                                    'in the new structure').format(i1,i2,block))
                        continue
                    if i1_sol[0]!=i2_sol[0]:
                        if show_warnings and np.max(np.abs(G[block][i1,i2].data)) > warning_threshold:
                            if mpi.is_master_node():
                                warn(('Element {},{} of block {} of G is approximated '+
                                    'to zero to match the new structure. Max abs value: {}').format(
                                        i1,i2,block,np.max(np.abs(G[block][i1,i2].data))))
                        continue
                    G_new[i1_sol[0]][i1_sol[1],i2_sol[1]] = \
                            G[block][i1,i2]
        return G_new

 def selfenergy(data, frozen_boson):
   if mpi.is_master_node():
     print "selfenergy: frozen_bozon: ",frozen_boson
   data.get_Sigma_loc_from_local_bubble()
   if not frozen_boson: data.get_P_loc_from_local_bubble()
   data.get_Sigmakw()
   data.get_Xkw() #if using optimized scheme make sure this is the order of calls (Sigmakw, Xkw then Pqnu)
   if not frozen_boson: data.get_Pqnu()
   if mpi.is_master_node():
     print "done with selfenergy"

 def monitor(self):
   try:  
     self.values.append(copy.deepcopy(self.mq()))    
   except:
     if mpi.is_master_node(): print "monitor: ",h5key," cound not read the value. appending nan..."
     self.values.append(float('nan'))    
   if mpi.is_master_node() and (not (self.archive_name is None)):
     A = HDFArchive(self.archive_name)
     A[self.h5key] = self.values
     del A
     print "monitor: ",self.h5key," : ", self.values[-1]

 def f(G1,G2) :
     #print abs(G1.data - G2.data)
     dS = max(abs(G1.data - G2.data).flatten())
     aS = max(abs(G1.data).flatten())
     if mpi.is_master_node():
         print "  Distances:", dS, " vs ", aS * dist
     return dS <= aS*dist

  def lattice(data, funcG, funcW, n): #the only option - we need Gkw and Wqnu for self energy in the next iteration
    #data.get_Gkw(funcG) #gets Gkw from G0 and Sigma
    def func(var, data):
      mu = var[0]
      dt = data[0]
      #print "func call! mu: ", mu, " n: ",dt.ns['up']
      n= data[1] 
      dt.mus['up'] = mu
      dt.mus['down'] = dt.mus['up']
      dt.get_Gkw_direct(funcG) #gets Gkw from w, mu, epsilon and Sigma and X
      dt.get_G_loc() #gets G_loc from Gkw
      dt.get_n_from_G_loc()     
      #print "funcvalue: ",-abs(n - dt.ns['up'])  
      return 1.0-abs(n - dt.ns['up'])  
    if not MASTER_SLAVE_ARCHITECTURE: mpi.barrier()
    varbest, funcvalue, iterations = amoeba(var=[data.mus['up']],
                                              scale=[0.01],
                                              func=func, 
                                              data = [data, n],
                                              itmax=30,
                                              ftolerance=1e-2,
                                              xtolerance=1e-2)
    if mpi.is_master_node():
      print "mu best: ", varbest
      print "-abs(diff n - data.n): ", funcvalue
      print "iterations used: ", iterations

    data.get_Gtildekw() #gets Gkw-G_loc

    data.get_Wqnu_from_func(funcW) #gets Wqnu from P and J 
    data.get_W_loc() #gets W_loc from Wqnu
    data.get_Wtildeqnu() #gets Wqnu-W_loc, 

 def __init__(self, archive, *args, **kwargs):
     self.archive = archive
     if not os.path.exists(archive) and mpi.is_master_node():
         archive = HDFArchive(archive, 'w')
         archive.create_group('results')
         archive['results']['n_dmft_loops'] = 0
         del archive
     mpi.barrier()

    def transport_coefficient(self, direction, iq, n, beta, method=None):
        r"""
        Calculates the transport coefficient A_n in a given direction for a given :math:`\Omega`. The required members (Gamma_w, directions, Om_mesh) have to be obtained first
        by calling the function :meth:`transport_distribution <pytriqs.applications.dft.sumk_dft_tools.SumkDFTTools.transport_distribution>`. For n>0 A is set to NaN if :math:`\Omega` is not 0.0. 
        
        Parameters
        ----------
        direction : string
           :math:`\alpha\beta` e.g.: 'xx','yy','zz','xy','xz','yz'.
        iq : integer
            Index of :math:`\Omega` point in the member Om_mesh.
        n : integer
            Number of the desired moment of the transport distribution.
        beta : double
            Inverse temperature :math:`\beta`.
        method : string
            Integration method: cubic spline and scipy.integrate.quad ('quad'), simpson rule ('simps'), trapezoidal rule ('trapz'), rectangular integration (otherwise)
            Note that the sampling points of the the self-energy are used!
        
        Returns
        -------
        A : double
            Transport coefficient.
        """

        if not (mpi.is_master_node()): return
        
        assert hasattr(self,'Gamma_w'), "transport_coefficient: Run transport_distribution first or load data from h5!"
        
        if (self.Om_mesh[iq] == 0.0 or  n == 0.0):
            A = 0.0
            # setup the integrand
            if (self.Om_mesh[iq] == 0.0):
                A_int = self.Gamma_w[direction][iq] * (self.fermi_dis(self.omega,beta) * self.fermi_dis(-self.omega,beta)) * (self.omega*beta)**n
            elif (n == 0.0):
                A_int = self.Gamma_w[direction][iq] * (self.fermi_dis(self.omega,beta) - self.fermi_dis(self.omega+self.Om_mesh[iq],beta))/(self.Om_mesh[iq]*beta)
            
            # w-integration
            if method == 'quad':
                # quad on interpolated w-points with cubic spline
                A_int_interp = interp1d(self.omega,A_int,kind='cubic')
                A = quad(A_int_interp, min(self.omega), max(self.omega), epsabs=1.0e-12,epsrel=1.0e-12,limit = 500)
                A = A[0]
            elif method == 'simps':
                # simpson rule for w-grid
                A = simps(A_int,self.omega)
            elif method == 'trapz':
                # trapezoidal rule for w-grid
                A = numpy.trapz(A_int,self.omega)
            else:
                # rectangular integration for w-grid (orignal implementation) 
                d_w = self.omega[1] - self.omega[0]
                for iw in xrange(self.Gamma_w[direction].shape[1]):
                    A += A_int[iw]*d_w
            A = A * numpy.pi * (2.0-self.SP)
        else:
            A = numpy.nan
        return A

    def conductivity_and_seebeck(self, beta, method=None):
        r"""
        Calculates the Seebeck coefficient and the optical conductivity by calling 
        :meth:`transport_coefficient <pytriqs.applications.dft.sumk_dft_tools.SumkDFTTools.transport_coefficient>`. 
        The required members (Gamma_w, directions, Om_mesh) have to be obtained first by calling the function 
        :meth:`transport_distribution <pytriqs.applications.dft.sumk_dft_tools.SumkDFTTools.transport_distribution>`. 

        Parameters
        ----------
        beta : double
            Inverse temperature :math:`\beta`.

        Returns
        -------
        optic_cond : dictionary of double vectors
            Optical conductivity in each direction and frequency given by Om_mesh.

        seebeck : dictionary of double
            Seebeck coefficient in each direction. If zero is not present in Om_mesh the Seebeck coefficient is set to NaN.
        """

        if not (mpi.is_master_node()):
            return

        assert hasattr(
            self, 'Gamma_w'), "conductivity_and_seebeck: Run transport_distribution first or load data from h5!"
        n_q = self.Gamma_w[self.directions[0]].shape[0]

        A0 = {direction: numpy.full((n_q,), numpy.nan)
              for direction in self.directions}
        A1 = {direction: numpy.full((n_q,), numpy.nan)
              for direction in self.directions}
        self.seebeck = {direction: numpy.nan for direction in self.directions}
        self.optic_cond = {direction: numpy.full(
            (n_q,), numpy.nan) for direction in self.directions}

        for direction in self.directions:
            for iq in xrange(n_q):
                A0[direction][iq] = self.transport_coefficient(
                    direction, iq=iq, n=0, beta=beta, method=method)
                A1[direction][iq] = self.transport_coefficient(
                    direction, iq=iq, n=1, beta=beta, method=method)
                print "A_0 in direction %s for Omega = %.2f    %e a.u." % (direction, self.Om_mesh[iq], A0[direction][iq])
                print "A_1 in direction %s for Omega = %.2f    %e a.u." % (direction, self.Om_mesh[iq], A1[direction][iq])
                if ~numpy.isnan(A1[direction][iq]):
                    # Seebeck is overwritten if there is more than one Omega =
                    # 0 in Om_mesh
                    self.seebeck[direction] = - \
                        A1[direction][iq] / A0[direction][iq] * 86.17
            self.optic_cond[direction] = beta * \
                A0[direction] * 10700.0 / numpy.pi
            for iq in xrange(n_q):
                print "Conductivity in direction %s for Omega = %.2f       %f  x 10^4 Ohm^-1 cm^-1" % (direction, self.Om_mesh[iq], self.optic_cond[direction][iq])
                if not (numpy.isnan(A1[direction][iq])):
                    print "Seebeck in direction      %s for Omega = 0.00      %f  x 10^(-6) V/K" % (direction, self.seebeck[direction])

        return self.optic_cond, self.seebeck

def spatial_FT_single_core(Qij):        
  if mpi.is_master_node(): print "spatial_FT_single core"
  n = len(Qij[:,0,0]) 
  nk = len(Qij[0,:,0])    
  Qk = numpy.zeros((n,nk,nk), dtype=numpy.complex_)
        
  Qk[:,:,:] = [ spf(Qij[l,:,:]) for l in range(n)]

  return Qk

 def __save_eal(self,Filename,it):
     if mpi.is_master_node():
         f=open(Filename,'a')
         f.write('\neff. atomic levels, Iteration %s\n'%it)
         for i in range(self.Nlm*self.Nspin):
             for j in range(self.Nlm*self.Nspin):
                 f.write("%10.6f %10.6f   "%(self.ealmat[i,j].real,self.ealmat[i,j].imag))
             f.write("\n")
         f.close()

 def __init__(self, **kwargs):
     self.parameters = dict()
     for i in kwargs:
         self.parameters[i] = kwargs[i]
     super(CDmft, self).__init__(**self.parameters)
     if mpi.is_master_node():
         archive = HDFArchive(self.parameters['archive'], 'a')
         archive['CDmft_version'] = CDmft._version
         del archive

    def save(self):
        """Saves some quantities into an HDF5 arxiv"""

        if not (mpi.is_master_node()): return # do nothing on nodes

        ar=HDFArchive(self.hdf_file,'a')
        ar[self.lda_data]['chemical_potential'] = self.chemical_potential
        ar[self.lda_data]['dc_energ'] = self.dc_energ
        ar[self.lda_data]['dc_imp'] = self.dc_imp
        del ar      

def spatial_FT(Qij, N_cores=1):   
  if N_cores == 1: return spatial_FT_single_core(Qij)     
  if mpi.is_master_node(): print "spatial_FT, N_cores: ",N_cores
  n = len(Qij[:,0,0]) 
  nk = len(Qij[0,:,0])    
  Qk = numpy.zeros((n,nk,nk), dtype=numpy.complex_)
        
  pool = Pool(processes=N_cores)              # start worker processes             
  Qk[:,:,:] = pool.map(spf, [ Qij[l,:,:] for l in range(n)])  
  pool.close()      
  return Qk

def is_vasp_running(vasp_pid):
    """
    Tests if VASP initial process is still alive.
    """
    pid_exists = False
    if mpi.is_master_node():
        try:
            os.kill(vasp_pid, 0)
        except OSError, e:
            pid_exists = e.errno == errno.EPERM
        else:
            pid_exists = True

 def Sigma( fermionic_struct, bosonic_struct, 
            Sigma, G, W, Lambda, 
            func,
            su2_symmetry, ising_decoupling, 
            p = {'0': 1.0, '1': 1.0, 'z': 1.0, 'c': 2.0},
            overwrite = True ):
   if mpi.is_master_node(): print "MSA: ",MASTER_SLAVE_ARCHITECTURE 
   for U in fermionic_struct.keys():
     if su2_symmetry and U!='up': continue      
     if overwrite: Sigma[U].fill(0.0) 
     for V in fermionic_struct.keys():            
       for A in bosonic_struct.keys():     
         if not (A in p.keys()):
           if mpi.is_master_node(): print "Sigma WARNING: skipping block ",A
           continue 
         if (U!=V and A!='+-')or((U==V)and(A=='+-')): continue
         m = -1.0
         if (A=='1' or A=='z') and (not ising_decoupling): m*=3.0
         #print "p[",A,"]: ", p[A], " m: ",m
         Sigma[U] += p[A] * m * ( func( G1 = partial(G, key=V),   G2 = partial(W, key=A),  Lambda = lambda wi1, wi2: Lambda(A, wi1, wi2)  ) )
   if su2_symmetry and ('down' in fermionic_struct.keys()): 
     Sigma['down'] = copy.deepcopy(Sigma['up'])

 def test_trilex(data):
   data.__class__ = trilex_data
   data.promote(data.n_iw/2, data.n_iw/2)
   solvers.cthyb.run(data, no_fermionic_bath=False, symmetrize_quantities=True, 
                           trilex=True, n_w_f=data.n_iw_f, n_w_b=data.n_iw_b,
                           n_cycles=20000, max_time=10*60, hartree_shift = 0.0 )
   data.get_chi3_imp()
   data.get_chi3tilde_imp()
   data.get_Lambda_imp()
   data.get_Sigma_test()
   data.get_P_test()
   if mpi.is_master_node():
     data.dump_test(suffix='-final')

      def func(var, data):
        mu = var[0]
        dt = data[0]
        #print "func call! mu: ", mu, " n: ",dt.ns['up']
        n= data[1] 
        dt.mus['up'] = mu
        if 'down' in dt.fermionic_struct.keys(): dt.mus['down'] = dt.mus['up']
        get_n(dt)        #print "funcvalue: ",-abs(n - dt.ns['up'])  

        val = 1.0-abs(n - dt.ns['up'])  
        if mpi.is_master_node(): print "amoeba func call: val = ",val
        if val != val: return -1e+6
        else: return val

  def __init__(self, monitored_quantity, accuracy=3e-5, func=None, struct=None, archive_name=None, h5key='diffs'):
    #monitored quantity needs to be a function returning an object in case the object is rewritten (changed address)
    self.mq = monitored_quantity

    self.accuracy = accuracy
    self.diffs = []
    self.func = func

    self.archive_name = archive_name
    self.struct = struct
    self.h5key = h5key

    if mpi.is_master_node(): print "converger initiialized: archive_name: %s h5key: %s accr: %s"%(archive_name,h5key,accuracy)