Support spin-orbit coupling

A spin-orbit flag is introduced in the constuctor.

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;

c++/pomerol_ed.hpp

@@ -71,7 +71,7 @@ namespace pomerol2triqs {
     std::set<Pomerol::ParticleIndex> computed_ops;
     std::unique_ptr<Pomerol::FieldOperatorContainer> ops_container;
 
-    Pomerol::Lattice init();
+    Pomerol::Lattice init(bool spin_orbit);
     Pomerol::ParticleIndex lookup_pomerol_index(indices_t const &i) const;
     std::set<Pomerol::ParticleIndex> 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);

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")
 

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

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)""",