Implement hubbard_one_solver for DMFT

Spin-orbit coupling not supported yet.

example/hubbard_one.py

@@ -0,0 +1,118 @@
+from pytriqs.applications.impurity_solvers.pomerol2triqs.hubbard_one_solver import Solver
+
+from pytriqs.archive import HDFArchive
+from pytriqs.gf 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
+from pytriqs.operators.util.hamiltonians import h_int_slater
+from pytriqs.operators.util.observables import N_op, S_op, L_op, LS_op
+from pytriqs.utility import mpi
+from itertools import product
+import numpy as np
+
+# Example of using hubbard_one_solver
+# 5-orbital atom with Slater interaction term
+
+####################
+# Input parameters #
+####################
+
+# Angular momentumd
+L = 2
+
+# Inverse temperature
+beta = 10.0
+
+# Slater parameters
+U = 5.0
+J = 0.1
+# F0 = U
+# F2 = J*(14.0/(1.0 + 0.625))
+# F4 = F2*0.625
+
+mu = U*0.5  # n=1
+# mu = U*1.5  # n=2
+
+# spin-orbit coupling
+ls_coupling = 0.0
+
+# Number of Matsubara frequencies for GF calculation
+n_iw = 200
+
+# Number of imaginary time slices for GF calculation
+n_tau = 1001
+
+# Number of bosonic Matsubara frequencies for G^2 calculations
+g2_n_wb = 5
+# Number of fermionic Matsubara frequencies for G^2 calculations
+g2_n_wf = 10
+
+#####################
+# Input for Pomerol #
+#####################
+
+spin_names = ("up", "dn")
+orb_names = range(-L, L+1)
+
+# GF structure
+off_diag = True
+gf_struct = set_operator_structure(spin_names, orb_names, off_diag=off_diag)
+
+# mkind = get_mkind(off_diag, None)
+# # Conversion from TRIQS to Pomerol notation for operator indices
+# index_converter = {mkind(sn, bn) : ("atom", bi, "down" if sn == "dn" else "up")
+#                    for sn, (bi, bn) in product(spin_names, enumerate(orb_names))}
+#
+# if mpi.is_master_node():
+#     print "Block structure of single-particle Green's functions:", gf_struct
+#     print "index_converter:", index_converter
+
+# Operators
+N = N_op(spin_names, orb_names, off_diag=off_diag)
+Sz = S_op('z', spin_names, orb_names, off_diag=off_diag)
+Lz = L_op('z', spin_names, orb_names, off_diag=off_diag, basis = 'spherical')
+Jz = Sz + Lz
+
+# LS coupling term
+# LS = ls_op(spin_names, orb_names, off_diag=off_diag, basis = 'spherical')
+
+# Hamiltonian
+# U_mat = U_matrix(L, radial_integrals = [F0,F2,F4], basis='spherical')
+U_mat = U_matrix(L, U_int=U, J_hund=J, basis='spherical')
+H_int = h_int_slater(spin_names, orb_names, U_mat, off_diag=off_diag)
+
+# H -= mu*N
+# H += ls_coupling * LS
+
+# Double check that we are actually using integrals of motion
+# h_comm = lambda op: H*op - op*H
+# str_if_commute = lambda op: "= 0" if h_comm(op).is_zero() else "!= 0"
+
+# print "Check integrals of motion:"
+# print "[H, N]", str_if_commute(N)
+# print "[H, Sz]", str_if_commute(Sz)
+# print "[H, Lz]", str_if_commute(Lz)
+# print "[H, Jz]", str_if_commute(Jz)
+
+E_levels = {}
+for sp in spin_names:
+    E_levels[sp] = np.diag( [-mu]*len(orb_names) )
+# print type(E_levels)
+# print type(E_levels['up'])
+
+const_of_motion = [N, Sz, Lz]
+
+#####################
+# Pomerol ED solver #
+#####################
+
+S = Solver(beta, gf_struct, 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")
+
+if mpi.is_master_node():
+    with HDFArchive('hubbard_one.out.h5', 'w') as ar:
+        ar['G_iw'] = S.G_iw
+        ar['G0_iw'] = S.G0_iw
+        ar['Sigma_iw'] = S.Sigma_iw

python/hubbard_one_solver.py

@@ -0,0 +1,97 @@
+
+from pomerol2triqs import PomerolED
+from pytriqs.gf.local import *
+from pytriqs.operators import Operator, c, c_dag, n
+from pytriqs.operators.util.op_struct import set_operator_structure, get_mkind
+import pytriqs.utility.mpi as mpi
+import numpy as np
+from itertools import product
+
+class Solver:
+    def __init__(self, beta, gf_struct, 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.verbose = verbose
+
+        if self.verbose and mpi.is_master_node():
+            print "\n*** Hubbard I solver using Pomerol library"
+            print "*** gf_struct =", gf_struct
+
+        # get spin_name and orb_names from gf_struct
+        self.__analyze_gf_struct(gf_struct)
+
+        if self.verbose and mpi.is_master_node():
+            print "*** spin_names =", self.spin_names
+            print "*** orb_names  =", self.orb_names
+
+        off_diag = True
+        mkind = get_mkind(off_diag, None)
+        # Conversion from TRIQS to Pomerol notation for operator indices
+        index_converter = {mkind(sn, bn) : ("atom", bi, "down" if sn == "dn" else "up")
+                           for sn, (bi, bn) in product(self.spin_names, enumerate(self.orb_names))}
+
+        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=""):
+
+        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
+
+        if const_of_motion is None:
+            self.__ed.diagonalize(H_0+H_int)
+        else:
+            self.__ed.diagonalize(H_0+H_int, const_of_motion)
+
+        # save data
+        if file_quantum_numbers:
+            self.__ed.save_quantum_numbers(file_quantum_numbers)
+        if file_eigenvalues:
+            self.__ed.save_eigenvalues(file_eigenvalues)
+
+        # set density-matrix cutoff
+        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)
+
+        # 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 -= E_list
+        self.Sigma_iw <<= self.G0_iw - inverse(self.G_iw)
+        self.G0_iw.invert()
+
+    def __copy_E_levels(self, E_levels):
+        """Copy E_levels after checking the data structure"""
+
+        # check data structure
+        assert( isinstance(E_levels, dict) )
+        for key, value in E_levels.items():
+            assert( isinstance(value, np.ndarray) )
+            assert( value.shape == (len(self.orb_names),len(self.orb_names)) )
+        # copy
+        self.E_levels = E_levels.copy()
+
+    def __analyze_gf_struct(self, gf_struct):
+        assert( isinstance(gf_struct, dict))
+        self.spin_names = gf_struct.keys()
+        self.orb_names = gf_struct.values()[0]
+        for value in gf_struct.values():
+            assert( self.orb_names == value )  # check if all spin blocks have the same orbitals