Init G_iw in constructor in hubbard_one_solver

j-otsuki · Mar 13, 2018 · 3faaf5e · 3faaf5e
1 parent abecc1e
commit 3faaf5e
Show file tree

Hide file tree

Showing 2 changed files with 17 additions and 14 deletions.
diff --git a/example/hubbard_one.py b/example/hubbard_one.py
@@ -1,7 +1,7 @@
 from pytriqs.applications.impurity_solvers.pomerol2triqs.hubbard_one_solver import Solver
 
 from pytriqs.archive import HDFArchive
-from pytriqs.gf import *
+from pytriqs.gf.local import *
 from pytriqs.operators import *
 from pytriqs.operators.util.op_struct import set_operator_structure, get_mkind
 from pytriqs.operators.util.U_matrix import U_matrix
@@ -107,9 +107,9 @@
 # Pomerol ED solver #
 #####################
 
-S = Solver(beta, gf_struct, spin_orbit=False, verbose=True)
+S = Solver(beta, gf_struct, n_iw, spin_orbit=False, verbose=True)
 
-S.solve(H_int, E_levels, n_iw, const_of_motion=const_of_motion, file_quantum_numbers="quantum_numbers.dat", file_eigenvalues="eigenvalues.dat")
+S.solve(H_int, E_levels, const_of_motion=const_of_motion, file_quantum_numbers="quantum_numbers.dat", file_eigenvalues="eigenvalues.dat")
 
 if mpi.is_master_node():
     with HDFArchive('hubbard_one.out.h5', 'w') as ar:

diff --git a/python/hubbard_one_solver.py b/python/hubbard_one_solver.py
@@ -8,19 +8,21 @@
 from itertools import product
 
 class Solver:
-    def __init__(self, beta, gf_struct, spin_orbit=False, verbose=True):
+    def __init__(self, beta, gf_struct, n_iw, spin_orbit=False, verbose=True):
 
         # TODO: spin_orbit
         if spin_orbit:
             print "*** hubbard_one_solver.Solver: spin_orbit not implemented yet"
 
         self.beta = beta
         self.gf_struct = gf_struct
+        self.n_iw = n_iw
         self.verbose = verbose
 
         if self.verbose and mpi.is_master_node():
             print "\n*** Hubbard I solver using Pomerol library"
             print "*** gf_struct =", gf_struct
+            print "*** n_iw =", n_iw
 
         # get spin_name and orb_names from gf_struct
         self.__analyze_gf_struct(gf_struct)
@@ -37,15 +39,20 @@ def __init__(self, beta, gf_struct, spin_orbit=False, verbose=True):
 
         self.__ed = PomerolED(index_converter, verbose)
 
-    def solve(self, H_int, E_levels, n_iw, const_of_motion=None, density_matrix_cutoff=1e-10, file_quantum_numbers="", file_eigenvalues=""):
+        # init G_iw
+        glist = lambda : [ GfImFreq(indices=self.orb_names, beta=beta, n_points=n_iw) for block, inner in gf_struct.items() ]
+        self.G_iw = BlockGf(name_list=self.spin_names, block_list=glist(), make_copies=False)
+        self.G_iw.zero()
+        self.G0_iw = self.G_iw.copy()
+        self.Sigma_iw = self.G_iw.copy()
+
+    def solve(self, H_int, E_levels, const_of_motion=None, density_matrix_cutoff=1e-10, file_quantum_numbers="", file_eigenvalues=""):
 
-        self.n_iw = n_iw
         self.__copy_E_levels(E_levels)
 
         H_0 = sum( self.E_levels[sn][o1,o2].real*c_dag(sn,o1)*c(sn,o2) for sn, o1, o2 in product(self.spin_names, self.orb_names, self.orb_names))
         if self.verbose and mpi.is_master_node():
             print "\n*** compute G_iw and Sigma_iw"
-            print "*** n_iw =", n_iw
             print "*** E_levels =", E_levels
             print "*** H_0 =", H_0
             print "*** const_of_motion =", const_of_motion
@@ -65,17 +72,13 @@ def solve(self, H_int, E_levels, n_iw, const_of_motion=None, density_matrix_cuto
         self.__ed.set_density_matrix_cutoff(density_matrix_cutoff)
 
         # Compute G(i\omega)
-        self.G_iw = self.__ed.G_iw(self.gf_struct, self.beta, self.n_iw)
-        # print type(self.G_iw)
+        self.G_iw << self.__ed.G_iw(self.gf_struct, self.beta, self.n_iw)
 
         # Compute G0 and Sigma
-        self.G0_iw = self.G_iw.copy()
-        self.Sigma_iw = self.G_iw.copy()
-
         E_list = [ self.E_levels[sn] for sn in self.spin_names ]  # from dict to list
-        self.G0_iw <<= iOmega_n
+        self.G0_iw << iOmega_n
         self.G0_iw -= E_list
-        self.Sigma_iw <<= self.G0_iw - inverse(self.G_iw)
+        self.Sigma_iw << self.G0_iw - inverse(self.G_iw)
         self.G0_iw.invert()
 
     def __copy_E_levels(self, E_levels):