From afc5ae0729d1046053092908a47f738cbea09334 Mon Sep 17 00:00:00 2001 From: Junya Otsuki Date: Fri, 9 Mar 2018 18:45:01 +0900 Subject: [PATCH] Implement hubbard_one_solver for DMFT Spin-orbit coupling not supported yet. --- example/hubbard_one.py | 118 +++++++++++++++++++++++++++++++++++ python/hubbard_one_solver.py | 97 ++++++++++++++++++++++++++++ 2 files changed, 215 insertions(+) create mode 100644 example/hubbard_one.py create mode 100644 python/hubbard_one_solver.py diff --git a/example/hubbard_one.py b/example/hubbard_one.py new file mode 100644 index 0000000..aa90a70 --- /dev/null +++ b/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 diff --git a/python/hubbard_one_solver.py b/python/hubbard_one_solver.py new file mode 100644 index 0000000..51b0a80 --- /dev/null +++ b/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