From 0eb1a21d35542cbc432bad9ec68e5518cceea12f Mon Sep 17 00:00:00 2001 From: Junya Otsuki Date: Mon, 19 Mar 2018 15:20:52 +0900 Subject: [PATCH] Support spin-orbit coupling A spin-orbit flag is introduced in the constuctor. --- c++/pomerol_ed.cpp | 9 +++--- c++/pomerol_ed.hpp | 4 +-- example/hubbard_one.py | 54 +++++++++++++++++++++--------------- python/hubbard_one_solver.py | 20 +++++++------ python/pomerol2triqs_desc.py | 2 +- 5 files changed, 50 insertions(+), 39 deletions(-) diff --git a/c++/pomerol_ed.cpp b/c++/pomerol_ed.cpp index a0180da..b89bccc 100644 --- a/c++/pomerol_ed.cpp +++ b/c++/pomerol_ed.cpp @@ -44,7 +44,7 @@ namespace pomerol2triqs { return neg_o ? std::conj(res) : res; } - Pomerol::Lattice pomerol_ed::init() { + Pomerol::Lattice pomerol_ed::init(bool spin_orbit) { Pomerol::Lattice l; @@ -61,13 +61,14 @@ namespace pomerol2triqs { it->second = std::max(it->second, pomerol_orb); } - for (auto const &site_orb : site_max_orb) l.addSite(new Pomerol::Lattice::Site(site_orb.first, site_orb.second + 1, 2)); + int n_spin = spin_orbit ? 1 : 2; + for (auto const &site_orb : site_max_orb) l.addSite(new Pomerol::Lattice::Site(site_orb.first, site_orb.second + 1, n_spin)); return l; } - pomerol_ed::pomerol_ed(index_converter_t const &index_converter, bool verbose) - : verbose(verbose), index_converter(index_converter), bare_lattice(init()), index_info(bare_lattice.getSiteMap()) { + pomerol_ed::pomerol_ed(index_converter_t const &index_converter, bool verbose, bool spin_orbit) + : verbose(verbose), index_converter(index_converter), bare_lattice(init(spin_orbit)), index_info(bare_lattice.getSiteMap()) { index_info.prepare(); if (verbose && !comm.rank()) { std::cout << "\nPomerol: lattice sites" << std::endl; diff --git a/c++/pomerol_ed.hpp b/c++/pomerol_ed.hpp index 8c10826..fa33505 100644 --- a/c++/pomerol_ed.hpp +++ b/c++/pomerol_ed.hpp @@ -71,7 +71,7 @@ namespace pomerol2triqs { std::set computed_ops; std::unique_ptr ops_container; - Pomerol::Lattice init(); + Pomerol::Lattice init(bool spin_orbit); Pomerol::ParticleIndex lookup_pomerol_index(indices_t const &i) const; std::set gf_struct_to_pomerol_indices(gf_struct_t const &gf_struct) const; double diagonalize_prepare(many_body_op_t const &hamiltonian); @@ -91,7 +91,7 @@ namespace pomerol2triqs { public: /// Create a new solver object - pomerol_ed(index_converter_t const &index_converter, bool verbose = false); + pomerol_ed(index_converter_t const &index_converter, bool verbose = false, bool spin_orbit = false); /// Diagonalize Hamiltonian optionally employing conservation of N and S_z void diagonalize(many_body_op_t const &hamiltonian, bool ignore_symmetries = false); diff --git a/example/hubbard_one.py b/example/hubbard_one.py index e9d3e5c..5af6b9a 100644 --- a/example/hubbard_one.py +++ b/example/hubbard_one.py @@ -52,26 +52,28 @@ # Input for Pomerol # ##################### -spin_names = ("up", "down") -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 +spin_names = ("up", "dn") +orb_names = map(str, range(-L, L+1)) +# orb_names = range(-L, L+1) # This also works + +flag_so = True if ls_coupling != 0 else False + +# GF and operator structure +if not flag_so: + # ('up', '-2'), ('dn', '-2'), etc + gf_struct = set_operator_structure(spin_names, orb_names, off_diag=True) + map_operator_structure = { (s,o) : (s,o) for s, o in product(spin_names, orb_names) } +else: + # ('ud', 'up_-2'), ('ud', 'dn_-2'), etc + gf_struct = {'ud': [str(s)+'_'+str(o) for s, o in product(spin_names, orb_names)] } + map_operator_structure = { (s,o) : ('ud', str(s)+'_'+str(o)) for s, o in product(spin_names, orb_names) } +# print gf_struct +# print map_operator_structure # 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') +N = N_op(spin_names, orb_names, map_operator_structure=map_operator_structure) +Sz = S_op('z', spin_names, orb_names, map_operator_structure=map_operator_structure) +Lz = L_op('z', spin_names, orb_names, map_operator_structure=map_operator_structure, basis = 'spherical') Jz = Sz + Lz # LS coupling term @@ -80,7 +82,7 @@ # 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_int = h_int_slater(spin_names, orb_names, U_mat, map_operator_structure=map_operator_structure) # H -= mu*N # H += ls_coupling * LS @@ -96,18 +98,24 @@ # print "[H, Jz]", str_if_commute(Jz) E_levels = {} -for sp in spin_names: - E_levels[sp] = np.diag( [-mu]*len(orb_names) ) +if not flag_so: + for sp in spin_names: + E_levels[sp] = np.diag( [-mu]*len(orb_names) ) +else: + # TODO: add matrix elements of spin-orbit coupling + E_levels['ud'] = np.diag( [-mu]*len(orb_names)*2 ) + # print type(E_levels) # print type(E_levels['up']) -const_of_motion = [N, Sz, Lz] +# const_of_motion = [N, Sz, Lz] +const_of_motion = [N, Sz, Lz] if not flag_so else [N, Jz] ##################### # Pomerol ED solver # ##################### -S = Solver(beta, gf_struct, n_iw, spin_orbit=False, verbose=True) +S = Solver(beta, gf_struct, n_iw, spin_orbit=flag_so, verbose=True) S.solve(H_int, E_levels, const_of_motion=const_of_motion, file_quantum_numbers="quantum_numbers.dat", file_eigenvalues="eigenvalues.dat") diff --git a/python/hubbard_one_solver.py b/python/hubbard_one_solver.py index d8e3819..8fb2f3d 100644 --- a/python/hubbard_one_solver.py +++ b/python/hubbard_one_solver.py @@ -10,10 +10,6 @@ class Solver: 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 @@ -32,13 +28,19 @@ def __init__(self, beta, gf_struct, n_iw, spin_orbit=False, verbose=True): print "*** spin_names =", self.spin_names print "*** orb_names =", self.orb_names - off_diag = True - mkind = get_mkind(off_diag, None) + # check spin_names + for sn in self.spin_names: + if not spin_orbit: assert sn in ('down', 'dn', 'up') # become either 'down' or 'up' + if spin_orbit: assert sn in ('down', 'dn', 'ud') # all become 'down' + # Conversion from TRIQS to Pomerol notation for operator indices - index_converter = {mkind(sn, bn) : ("atom", bi, "down" if sn == "dn" else sn) + # TRIQS: ('dn', '-2') --> Pomerol: ('atom', 0, 'down') + # NOTE: When spin_orbit is true, only spin 'down' is used. + mkind = get_mkind(True, None) + index_converter = {mkind(sn, bn) : ("atom", bi, "down" if sn == "dn" or sn == "ud" else sn) for sn, (bi, bn) in product(self.spin_names, enumerate(self.orb_names))} - self.__ed = PomerolED(index_converter, verbose) + self.__ed = PomerolED(index_converter, verbose, spin_orbit) # init G_iw glist = lambda : [ GfImFreq(indices=self.orb_names, beta=beta, n_points=n_iw) for block, inner in gf_struct.items() ] @@ -51,7 +53,7 @@ def solve(self, H_int, E_levels, const_of_motion=None, density_matrix_cutoff=1e- 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)) + H_0 = sum( self.E_levels[sn][oi1,oi2].real*c_dag(sn,on1)*c(sn,on2) for sn, (oi1,on1), (oi2,on2) in product(self.spin_names, enumerate(self.orb_names), enumerate(self.orb_names))) if self.verbose and mpi.is_master_node(): print "\n*** compute G_iw and Sigma_iw" print "*** E_levels =", E_levels diff --git a/python/pomerol2triqs_desc.py b/python/pomerol2triqs_desc.py index 8c99de5..0e41211 100644 --- a/python/pomerol2triqs_desc.py +++ b/python/pomerol2triqs_desc.py @@ -40,7 +40,7 @@ doc = r"", # doc of the C++ class ) -c.add_constructor("""(index_converter_t index_converter, bool verbose = false)""", +c.add_constructor("""(index_converter_t index_converter, bool verbose = false, bool spin_orbit = false)""", doc = """Create a new solver object """) c.add_method("""void diagonalize (many_body_op_t hamiltonian, bool ignore_symmetries = false)""",