From 3b60e530199ddd456298dd6025b2cd4824b7f03f Mon Sep 17 00:00:00 2001
From: Dennis Heffels <d.heffels@fz-juelich.de>
Date: Wed, 20 Jan 2021 10:21:27 +0000
Subject: [PATCH] phi_matrix

---
 kwant/builder.py       | 28 +++++++++++++++++++++-------
 kwant/physics/gauge.py | 23 ++++++++++-------------
 2 files changed, 31 insertions(+), 20 deletions(-)

diff --git a/kwant/builder.py b/kwant/builder.py
index e01980bd..f65a8e08 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -1528,13 +1528,24 @@ class Builder:
     _require_system
 
 
-def _add_peierls_phase(syst, q, peierls_parameter):
+def _add_peierls_phase(syst, peierls_parameter, signs):
+    def phi_matrix(phi, signs):
+        phase_factors = []
+        for sign in signs:
+            if sign == -1:
+                phase_factors.append(phi.conjugate())
+            elif sign == 1:
+                phase_factors.append(phi)
+        return np.diag(phase_factors)
 
     @memoize
     def wrap_hopping(hop):
         def f(*args):
             a, b, *args, phi = args
-            return hop(a, b, *args) * phi(a, b)
+            if signs is not None:
+                return np.dot(hop(a, b, *args), phi_matrix(phi(a, b), signs))
+            else:
+                return hop(a, b, *args) * phi(a, b)
 
         params = ('_site0', '_site1')
         if callable(hop):
@@ -1550,7 +1561,10 @@ def _add_peierls_phase(syst, q, peierls_parameter):
     const_hoppings = {}
 
     def const_hopping(a, b, phi):
-        return const_hoppings[a, b] * phi(a, b)
+        if signs is not None:
+            return np.dot(hop(a, b, *args), phi_matrix(phi(a, b), signs))
+        else:
+            return hop(a, b, *args) * phi(a, b)
 
     # fix the signature
     params = ('_site0', '_site1', peierls_parameter)
@@ -1570,8 +1584,8 @@ def _add_peierls_phase(syst, q, peierls_parameter):
             BuilderLead(
                 _add_peierls_phase(
                     lead.builder,
-                    q,
                     peierls_parameter + '_lead{}'.format(i),
+                    signs
                 ),
                 lead.interface,
                 lead.padding,
@@ -1591,7 +1605,7 @@ def _add_peierls_phase(syst, q, peierls_parameter):
     return ret
 
 
-def add_peierls_phase(syst, q=1, peierls_parameter='phi', fix_gauge=True):
+def add_peierls_phase(syst, peierls_parameter='phi', fix_gauge=True, signs=None):
     """Add a Peierls phase parameter to a Builder.
 
     Parameters
@@ -1691,10 +1705,10 @@ def add_peierls_phase(syst, q=1, peierls_parameter='phi', fix_gauge=True):
         raise TypeError("'add_peierls_phase' does not work with "
                         "vectorized Builders")
 
-    ret = _add_peierls_phase(syst, q, peierls_parameter).finalized()
+    ret = _add_peierls_phase(syst, peierls_parameter, signs).finalized()
 
     if fix_gauge:
-        return ret, wrap_gauge(magnetic_gauge(ret, q))
+        return ret, wrap_gauge(magnetic_gauge(ret))
     else:
         return ret
 
diff --git a/kwant/physics/gauge.py b/kwant/physics/gauge.py
index 1ffee248..286b75f3 100644
--- a/kwant/physics/gauge.py
+++ b/kwant/physics/gauge.py
@@ -760,7 +760,7 @@ def _check_composite_system(syst):
 
 ### Phase calculation
 
-def _calculate_phases(loops, pos, previous_phase, flux, q):
+def _calculate_phases(loops, pos, previous_phase, flux):
     """Calculate the phase across the terminal links of a set of loops
 
     Parameters
@@ -788,7 +788,7 @@ def _calculate_phases(loops, pos, previous_phase, flux, q):
     for loop in loops:
         tail, head = loop[-1], loop[0]
         integral = flux([pos(p) for p in loop])
-        phase = np.exp(q* 1j * np.pi * integral)
+        phase = np.exp(1j * np.pi * integral)
         phases[tail, head] = phase / previous_phase(phases, loop)
     return phases
 
@@ -879,7 +879,7 @@ def _infinite_wrapper(syst, phases, a, b):
         return 1
 
 
-def _peierls_finite(syst, loops, q, syst_field, tol, average):
+def _peierls_finite(syst, loops, syst_field, tol, average):
     integrate = partial(_surface_integral, syst_field,
                         tol=tol, average=average)
     phases = _calculate_phases(
@@ -887,12 +887,11 @@ def _peierls_finite(syst, loops, q, syst_field, tol, average):
         syst.pos,
         _previous_phase_finite,
         integrate,
-        q,
     )
     return partial(_finite_wrapper, syst, phases)
 
 
-def _peierls_infinite(syst, loops, q, extended_sites, syst_field, tol, average):
+def _peierls_infinite(syst, loops, extended_sites, syst_field, tol, average):
     integrate = partial(_surface_integral, syst_field,
                         tol=tol, average=average)
     phases = _calculate_phases(
@@ -900,12 +899,11 @@ def _peierls_infinite(syst, loops, q, extended_sites, syst_field, tol, average):
         lambda i: extended_sites(i).pos,
         partial(_previous_phase_infinite, syst.cell_size),
         integrate,
-        q,
     )
     return partial(_infinite_wrapper, syst, phases)
 
 
-def _peierls_composite(syst, loops, q, which_patch, extended_sites, lead_gauges,
+def _peierls_composite(syst, loops, which_patch, extended_sites, lead_gauges,
                        syst_field, *lead_fields, tol, average):
     if len(lead_fields) != len(syst.leads):
         raise ValueError('Magnetic fields must be provided for all leads.')
@@ -925,8 +923,7 @@ def _peierls_composite(syst, loops, q, which_patch, extended_sites, lead_gauges,
         lambda i: extended_sites(i).pos,
         partial(_previous_phase_composite,
                 which_patch, extended_sites, lead_phases),
-        flux,
-        q,
+        flux
     )
 
     return (partial(_finite_wrapper, syst, phases), *lead_phases)
@@ -977,21 +974,21 @@ class magnetic_gauge:
     >>> kwant.hamiltonian_submatrix(syst, params=params)
     """
 
-    def __init__(self, syst, q=1):
+    def __init__(self, syst):
         if isinstance(syst, builder.FiniteSystem):
             if syst.leads:
                 loops, which_patch, extended_sites = _loops_in_composite(syst)
                 lead_gauges = [magnetic_gauge(lead) for lead in syst.leads]
                 self._peierls = partial(_peierls_composite, syst,
-                                        loops, q, which_patch,
+                                        loops, which_patch,
                                         extended_sites, lead_gauges)
             else:
                 loops = _loops_in_finite(syst)
-                self._peierls = partial(_peierls_finite, syst, loops, q)
+                self._peierls = partial(_peierls_finite, syst, loops)
         elif isinstance(syst, builder.InfiniteSystem):
             loops, extended_sites = _loops_in_infinite(syst)
             self._peierls = partial(_peierls_infinite, syst,
-                                    loops, q, extended_sites)
+                                    loops, extended_sites)
         else:
             raise TypeError('Expected a finalized Builder')
 
-- 
GitLab