TRIQS / triqs

a Toolbox for Research on Interacting Quantum Systems
https://triqs.github.io
GNU General Public License v3.0
134 stars 69 forks source link

serial Hilbert transform #376

Open JaksaVucicevic opened 7 years ago

JaksaVucicevic commented 7 years ago

The Hilbert_transform function does not have an option for serial execution. I encountered a problem with this because I'm doing it in only one of the mpi processes, while the others are blocked by mpi.recv. What happens then is that Hilbert transform waits forever at reduce.

I propose the following solution at pytriqs/dos/hilbert_transform.py:69:

 def __call__ (self, Sigma, mu=0, eta=0, field=None, epsilon_hat=None, result=None,
                  n_points_integral=None, test_convergence=None, mpi_parallel = True):
        #...not showing all the code
        def HT(res):
            #... not showing all the code
            dehe = [self.rho_for_sum, eps_hat, self.dos.eps]
            for d, e_h, e in ( (itertools.izip (*dehe)) if (not mpi_parallel) else itertools.izip(*[mpi.slice_array(A) for A in dehe]) ):
                tmp2.copy_from(tmp)
                tmp2 -= e_h
                if Sigma_fnt: tmp2 -= Sigma(e)
                tmp2.invert()
                tmp2 *= d
                res += tmp2
            # sum the res GF of all nodes and returns the results on all nodes...
            # Cf Boost.mpi.python, collective communicator for documentation.
            # The point is that res is pickable, hence can be transmitted between nodes without further code...
            if mpi_parallel: 
                res << mpi.all_reduce(mpi.world, res, lambda x, y: x+y)
                mpi.barrier()

         #...... rest of the code where HT is being called

I checked that it works.

Wentzell commented 7 years ago

Can you post a minimal example to reproduce the error?

JaksaVucicevic commented 7 years ago

The following works in serial execution, but freezes after it's done if you run it with

$ mpirun -np 2 pytriqs script.py

script.py:

from pytriqs.lattice.tight_binding import *
from pytriqs.dos.hilbert_transform import *
import pytriqs.utility.mpi as mpi

t    = -1.00                
U    = 2
beta = 10
n_iw = 100
mu = U/2.0 + 0.1

BL = BravaisLattice(units = [(1,0,0) , (0,1,0) ])

hop= {  (1,0)  :  [[ t]],       
        (-1,0) :  [[ t]],     
        (0,1)  :  [[ t]],
        (0,-1) :  [[ t]]}

TB = TightBinding ( BL, hop)

ht = HilbertTransform(dos(TB, n_kpts= 500, n_eps = 101, name = 'dos')[0])

G_iw = GfImFreq(indices = [0], beta = beta, n_points = n_iw, statistic = 'Fermion')
Sigma_iw = G_iw.copy()
Sigma_iw << U/2.0

if mpi.is_master_node():
    print '[ Node %s]:'%mpi.rank,'doing the Hilbert transform'
    G_iw = ht(Sigma = Sigma_iw, mu = mu)
else:
    print '[ Node %s]:'%mpi.rank,'doing something else'
pram2010 commented 3 years ago

Hi, I am using the HilbertTransform (#655) from the triqs library. Bellow is some segment of the code.

Argument in HilbertTransform >> epsilon_hat

Without putting this argument, I am getting the correct result. But, I want set the value of epsilon_hat for my calculation. It will be great if you can assist me. Thanks.

Code:

from pytriqs.lattice import BravaisLattice, BrillouinZone from pytriqs.gf import from pytriqs.dos import

import warnings warnings.filterwarnings("ignore")

import numpy as np from math import cos, pi ########################################

def dos_circ(eps): return 2.0/(np.pi)*np.sqrt(1.0-eps**2)

x_max = 0.9999 x_min = - x_max d = DOSFromFunction(dos_circ, x_min, x_max)

H = HilbertTransform(d)

e_max=3 G = GfReFreq(indices = [0,1], window=(-e_max,e_max)) Sigma0 = G.copy()

G << H(Sigma = Sigma0, mu=0., eta=0.05) # working

G << H(Sigma = Sigma0, mu=0., eta=0.05, epsilon_hat=??)