major rewrite: fully include AC voltage in RG flow

README.md

@@ -34,9 +34,7 @@ The experimental code for `rtrg_cublas` aims to provide the same functionality a
 Note that `rtrg_cublas` is probably less stable and can be much slower than `rtrg_c` if not configured correctly.
 
 ### Floquet matrix classes
-The classes `RGfunction` and `SymRGfunction` in `rtrg.py` and `compact_rtrg.py` provide a high-level interface for calculations with Floquet matrices in python.
-Both were written specifically for the Floquet real-time renormalization group (FRTRG), but the basis of these classes may be more generally applicable.
-`SymRGfunction` is written for the special case of Floquet matrices representing functions in time-domain which fulfil f(t+T/2) = ±f(t) where T is the period of the system.
+The class `RGfunction` in `rtrg.py` provide a high-level interface for calculations with Floquet matrices including replicas with shifted energy argument in python.
 
 `reservoirmatrix.py` provides a matrix of `RGfunction` objects and some functions that speed up some calculations with these objects by using symmetries.
 

reservoirmatrix.py

@@ -13,6 +13,10 @@ class ReservoirMatrix:
     2x2 matrix of RGfunctions.
     This includes a system of book keeping for energy shifts by multiples of
     the voltage in products of ReservoirMatrices.
+    if symmetry != 0:
+        self.data[1,0] == symmetry * self.data[0,1].floquetConjugate()
+    if symmety != 0 and global_properties.symmetric:
+        self.data[0,0] == self.data[1,1]
     '''
 
     def __init__(self, global_properties, symmetry=0):
@@ -143,7 +147,7 @@ class ReservoirMatrix:
             res_01_10 = self[0,1] @ other[1,0]
             res[0,0] = res_00_00 + res_01_10
             res[0,1] = self[0,0] @ other[0,1] + self[0,1] @ other[1,1]
-            symmetry = self.symmetry * other.symmetry
+            symmetry = self.symmetry * other.symmetry * (self.voltage_shifts == 0)
             if symmetry == 0 or IGNORE_SYMMETRIES:
                 res[1,1] = self[1,0] @ other[0,1] + self[1,1] @ other[1,1]
                 res[1,0] = self[1,0] @ other[0,0] + self[1,1] @ other[1,0]
@@ -159,7 +163,7 @@ class ReservoirMatrix:
             # 3 multiplications with symmetry but xL != xR
             # 2 multiplications with symmetry and xL == xR
             assert self.global_properties is other.global_properties
-            res = ReservoirMatrix(self.global_properties, self.symmetry * other.symmetry)
+            res = ReservoirMatrix(self.global_properties, self.symmetry * other.symmetry * (self.voltage_shifts == 0))
             res.voltage_shifts = self.voltage_shifts + other.voltage_shifts
             res[0,0] = self[0,0] @ other
             if res.symmetry == 0 or IGNORE_SYMMETRIES:
@@ -181,7 +185,7 @@ class ReservoirMatrix:
         if isinstance(other, RGfunction):
             assert self.voltage_shifts == 0
             assert self.global_properties is other.global_properties
-            res = ReservoirMatrix(self.global_properties, self.symmetry * other.symmetry)
+            res = ReservoirMatrix(self.global_properties, self.symmetry * other.symmetry * (other.voltage_shifts == 0))
             res.voltage_shifts = other.voltage_shifts
             res[0,0] = other @ self[0,0]
             if res.symmetry == 0 or IGNORE_SYMMETRIES:
@@ -261,15 +265,22 @@ class ReservoirMatrix:
         return shifted
 
     def to_numpy_array(self):
-        array = np.ndarray((2,2,*self.shape()), dtype=np.complex128)
+        array = np.ndarray((2,2,*self[0,0].values.shape), dtype=np.complex128)
         for i in range(2):
             for j in range(2):
-                if isinstance(self[i,j], SymRGfunction):
-                    array[i,j] = self[i,j].toRGfunction().values
-                else:
                 array[i,j] = self[i,j].values
         return array
 
+    def check_symmetry(self):
+        assert self.symmetry in (-1,0,1)
+        if self.global_properties.symmetric:
+            assert np.allclose(self[0,0].values, self[1,1].values)
+        if self.symmetry:
+            assert np.allclose(self[0,1].values, self.symmetry*self[1,0].floquetConjugate().values)
+            self[0,0].check_symmetry()
+            self[1,1].check_symmetry()
+
+
 
 def einsum_34_12_43(a, b, c):
     '''

rtrg.py

@@ -36,9 +36,8 @@ class GlobalRGproperties:
             nmax = 0,
             padding = 0,
             voltage_branches = None,
-            resonant_dc_shift = 0,
             symmetric = True,
-            vdc = 0,
+            mu = 0,
             clear_corners = 0,
             )
 
@@ -49,7 +48,6 @@ class GlobalRGproperties:
         omega
         nmax
         voltage_branches = None
-        vdc = 0
         '''
         #print('Created new GlobalRGproperties object', flush=True)
         self.__dict__.update(kwargs)
@@ -71,11 +69,8 @@ class GlobalRGproperties:
         '''
         floquet_dim = 2*self.__dict__['nmax'] + 1
         try:
-            return (*self.__dict__['energies'].shape[:-1], floquet_dim, floquet_dim)
+            return (*self.__dict__['energies'].shape[:-1], 2*self.__dict__['voltage_branches']+1, floquet_dim, floquet_dim)
         except KeyError:
-            if self.__dict__['voltage_branches'] is None:
-                return (floquet_dim, floquet_dim)
-            else:
             return (2*self.__dict__['voltage_branches'] + 1, floquet_dim, floquet_dim)
 
     def copy(self):
@@ -110,8 +105,8 @@ class RGfunction:
 
     self.values:
     This contains a arrays of the shapes
-    (..., 2*nmax+1) as energies and
-    (..., 2*nmax+1, 2*nmax+1) as values.
+    (2*voltage_branches+1, 2*nmax+1) as energies and
+    (2*voltage_branches+1, 2*nmax+1, 2*nmax+1) as values.
 
     self.voltage_shifts:
     energy shifts (in units of DC voltage) which should be included in all
@@ -125,7 +120,7 @@ class RGfunction:
 
     self.symmetry:
     Valid values are 0, -1, +1. If non-zero, it states the symmetry
-    X[::-1,::-1] == symmetry * X.conjugate() for this Floquet matrix.
+    X[::-1,::-1,::-1] == symmetry * X.conjugate() for this Floquet matrix.
     '''
     def __init__(self, global_properties, values=None, voltage_shifts=0, symmetry=0, **kwargs):
         self.global_properties = global_properties
@@ -139,14 +134,10 @@ class RGfunction:
             self.values = self.initial()
         elif (type(values) == str and values == 'identity'):
             self.symmetry = 1
-            if self.global_properties.energies.ndim == 1:
-                self.values = np.identity(2*self.nmax+1, dtype=np.complex128).T
-            else:
             self.values = np.broadcast_to( np.identity(2*self.nmax+1, dtype=np.complex128).reshape(1, 2*self.nmax+1, 2*self.nmax+1), self.shape())
         else:
             self.values = np.asarray(values, dtype=np.complex128)
             assert self.values.shape[-1] == self.values.shape[-2] == 2*self.nmax+1 == self.energies.shape[-1]
-        assert self.values.ndim == self.energies.ndim + 1
 
     def __getattr__(self, name):
         return getattr(self.global_properties, name)
@@ -185,7 +176,7 @@ class RGfunction:
             assert abs(self.energies[self.nmax].real) < 1e-12
             return RGfunction(self.global_properties, np.conjugate(self.values[::-1,::-1]), self.voltage_shifts)
         elif self.values.ndim == 3:
-            assert abs(self.energies[self.voltage_branches, self.nmax].real) < 1e-12
+            assert abs(self.energies[self.nmax].real) < 1e-12
             return RGfunction(self.global_properties, np.conjugate(self.values[::-1,::-1,::-1]), self.voltage_shifts)
         else:
             raise NotImplementedError()
@@ -198,14 +189,17 @@ class RGfunction:
         Note: This is only approximately associative, as long the function
         converges to 0 for |n| of order of nmax.
         '''
-        if isinstance(other, SymRGfunction):
-            return NotImplemented
         if not isinstance(other, RGfunction):
             return NotImplemented
-        assert self.values.shape == other.values.shape
         assert self.global_properties is other.global_properties
+        if self.values.ndim == 2 and other.values.ndim == 3:
+            symmetry = self.symmetry * other.symmetry * (self.voltage_shifts == 0)
+            vb = other.values.shape[0]//2
+            matrix = rtrg_c.multiply_extended(other.values[vb+self.voltage_shifts].T, self.values.T, self.padding, symmetry, self.clear_corners).T
+        else:
+            assert self.values.shape == other.values.shape
             right = other.shift_energies(self.voltage_shifts)
-        symmetry = self.symmetry * other.symmetry * (self.voltage_shifts == 0) * (other.voltage_shifts == 0)
+            symmetry = self.symmetry * other.symmetry * (self.voltage_shifts == 0)
             matrix = rtrg_c.multiply_extended(right.values.T, self.values.T, self.padding, symmetry*(self.values.ndim==2), self.clear_corners).T
         if VERBOSE:
             RGfunction.MM_COUNTER[(symmetry != 0) | (self.values.T.flags.c_contiguous << 3) | (self.values.T.flags.f_contiguous << 4) | (right.values.T.flags.c_contiguous << 1) | (right.values.T.flags.f_contiguous << 2)]  += 1
@@ -213,13 +207,16 @@ class RGfunction:
         return res
 
     def __imatmul__(self, other):
-        if isinstance(other, SymRGfunction):
-            return NotImplemented
         if not isinstance(other, RGfunction):
             return NotImplemented
-        assert self.values.shape == other.values.shape
         assert self.global_properties is other.global_properties
         self.energy_shifted_copies.clear()
+        if self.values.ndim == 2 and other.values.ndim == 3:
+            symmetry = self.symmetry * other.symmetry * (self.voltage_shifts == 0)
+            vb = other.values.shape[0]//2
+            matrix = rtrg_c.multiply_extended(other.values[vb+self.voltage_shifts].T, self.values.T, self.padding, symmetry, self.clear_corners).T
+        else:
+            assert self.values.shape == other.values.shape
             right = other.shift_energies(self.voltage_shifts)
             self.symmetry *= other.symmetry
             self.values = rtrg_c.multiply_extended(right.values.T, self.values.T, self.padding, self.symmetry*(self.values.ndim==2), self.clear_corners, OVERWRITE_LEFT).T
@@ -411,8 +408,9 @@ class RGfunction:
     def __eq__(self, other):
         return ( self.global_properties is other.global_properties ) and self.voltage_shifts == other.voltage_shifts and np.allclose(self.values, other.values) and self.symmetry == other.symmetry
 
-    def k2lambda(self):
+    def k2lambda(self, shift_matrix=0):
         '''
+        shift is usually n*zinv*mu.
         Assume that self is K_n^m(E) = K_n(E-mΩ).
         Then calculate Λ_n^m(E) such that (approximately)
 
@@ -425,10 +423,50 @@ class RGfunction:
         overdetermined. This means that we can in general only get an
         approximate solution because an exact solution does not exist.
         '''
-        invert = -self
-        invert += self.energies
-        invert.symmetry = -1 if self.symmetry == -1 else 0
-        return invert.inverse()
+        assert self.values.ndim == 3
+        assert self.energies.ndim == 1
+        invert_array = np.ndarray((2*self.voltage_branches+5,2*self.nmax+1,2*self.nmax+1), dtype=np.complex128)
+        invert_array[2:-2] = -self.values + np.diag(self.energies).reshape((1,2*self.nmax+1,2*self.nmax+1)) + shift_matrix.values
+        symmetry = -(self.symmetry == -1)*(shift_matrix.symmetry == -1)
+        dim0 = 2*self.voltage_branches+1
+        if EXTRAPOLATE_VOLTAGE:
+            try:
+                interp = interp1d(np.arange(dim0), invert_array[2:-2], 'quadratic', axis=0, fill_value='extrapolate')
+            except ValueError:
+                interp = interp1d(np.arange(dim0), invert_array[2:-2], 'linear', axis=0, fill_value='extrapolate')
+            invert_array[-2:] = interp(dim0 + np.arange(2))
+            invert_array[:2] = interp(np.arange(-2, 0))
+        else:
+            invert_array[-2:] = invert_array[-3]
+            invert_array[-2:] += shift_matrix[shift_matrix.values.shape[0]//2+1:shift_matrix.values.shape[0]//2+3]
+            invert_array[:2] = invert_array[2]
+            invert_array[:2] += shift_matrix[shift_matrix.values.shape[0]//2-2:shift_matrix.values.shape[0]//2]
+        inverted = rtrg_c.invert_extended(invert_array.T, self.padding, round(LAZY_INVERSE_FACTOR*self.padding)).T
+        if VERBOSE:
+            RGfunction.INV_COUNTER[invert_array.T.flags.c_contiguous | (invert_array.T.flags.f_contiguous << 1)] += 1
+
+        res = RGfunction(self.global_properties, inverted[2:-2], voltage_shifts=self.voltage_shifts, symmetry=symmetry)
+        res.energy_shifted_copies[2] = RGfunction(self.global_properties, inverted[4:], self.voltage_shifts, symmetry=0)
+        res.energy_shifted_copies[1] = RGfunction(self.global_properties, inverted[3:-1], self.voltage_shifts, symmetry=0)
+        res.energy_shifted_copies[-1] = RGfunction(self.global_properties, inverted[1:-3], self.voltage_shifts, symmetry=0)
+        res.energy_shifted_copies[-2] = RGfunction(self.global_properties, inverted[:-4], self.voltage_shifts, symmetry=0)
+        return res
+
+    def reduced(self, shift=0):
+        '''
+        Remove voltage-shifted copies
+        '''
+        assert self.values.ndim == 3
+        assert abs(shift) <= self.voltage_branches
+        return RGfunction(self.global_properties, self.values[self.voltage_branches+shift], symmetry=self.symmetry*(shift==0))
+
+    def reduced_to_voltage_branches(self, voltage_branches):
+        if 2*voltage_branches+1 == self.values.shape[0]:
+            return self
+        assert 0 < voltage_branches < self.values.shape[0]//2
+        assert self.values.ndim == 3
+        diff = self.values.shape[0]//2 - voltage_branches
+        return RGfunction(self.global_properties, self.values[diff:-diff], symmetry=self.symmetry, voltage_shifts=self.voltage_shifts)
 
     def inverse(self):
         '''
@@ -442,6 +480,9 @@ class RGfunction:
         '''
 
         assert self.voltage_shifts == 0
+        if self.padding == 0:
+            res = np.linalg.inv(self.values)
+        else:
             try:
                 res = rtrg_c.invert_extended(self.values.T, self.padding, round(LAZY_INVERSE_FACTOR*self.padding)).T
                 if VERBOSE:
@@ -453,7 +494,6 @@ class RGfunction:
         return RGfunction(self.global_properties, res, self.voltage_shifts, self.symmetry)
 
     def shift_energies(self, n=0):
-        # TODO: use symmetries
         '''
         Shift energies by n*self.vdc. The calculated RGfunction is kept in cache.
         Assumptions:
@@ -462,24 +502,14 @@ class RGfunction:
         * derivative is a RGfunction or None
         * if the length of the last axis of self.energies is < 2, then derivative must not be None
         '''
-        if n==0 or self.vdc==0:
+        if n==0:
             return self
         try:
             return self.energy_shifted_copies[n]
         except KeyError:
             assert self.values.ndim == 3
             newvalues = np.ndarray(self.values.shape, dtype=np.complex128)
-            try:
-                assert self.derivative.global_properties is self.global_properties
-                if n > 0:
-                    newvalues[:-n] = self.values[n:]
-                    for i in range(n):
-                        newvalues[i-n] = self.values[-1] + i*self.vdc*self.derivative[-1]
-                else:
-                    newvalues[-n:] = self.values[:n]
-                    for i in range(-n):
-                        newvalues[i] = self.values[0] - (n-i)*self.vdc*self.derivative[-1]
-            except AttributeError:
+            if EXTRAPOLATE_VOLTAGE == 1:
                 try:
                     interp = interp1d(np.arange(self.values.shape[0]), self.values, 'quadratic', axis=0, fill_value='extrapolate')
                 except ValueError:
@@ -490,7 +520,21 @@ class RGfunction:
                 else:
                     newvalues[-n:] = self.values[:n]
                     newvalues[:-n] = interp(np.arange(n, 0))
+            else:
+                if n > 0:
+                    newvalues[:-n] = self.values[n:]
+                    newvalues[-n:] = self.values[-1:]
+                else:
+                    newvalues[-n:] = self.values[:n]
+                    newvalues[:-n] = self.values[:1]
             self.energy_shifted_copies[n] = RGfunction(self.global_properties, newvalues, self.voltage_shifts)
             return self.energy_shifted_copies[n]
 
-from compact_rtrg import SymRGfunction, OVERWRITE_LEFT, OVERWRITE_BOTH, OVERWRITE_RIGHT
+    def check_symmetry(self):
+        assert self.symmetry in (-1,0,1)
+        if self.symmetry:
+            if self.values.ndim == 2:
+                conjugate = self.values[::-1,::-1].conjugate()
+            elif self.values.ndim == 3:
+                conjugate = self.values[::-1,::-1,::-1].conjugate()
+            assert np.allclose(self.values, self.symmetry*conjugate)

settings.py

@@ -16,15 +16,27 @@ defaults = dict(
     # Log progress to stdout every LOG_TIME seconds
     LOG_TIME = 10, # in s
 
+    # Number of parallel threads (python multiprocessing)
+    PY_THREADS = 0,
+
     # Improve accuracy of current and conductance by terms which vanish in the limit D→∞
     EXTRA_INIT_CONDITIONS = 1,
 
     # Raise exception if no symmetries can be used in calculation steps.
     ENFORCE_SYMMETRIC = 0,
 
+    # Check symmetries before each iteration of the RG equations.
+    CHECK_SYMMETRIES = 0,
+
     # Do not use any symmetries.
     IGNORE_SYMMETRIES = 0,
 
+    # How to extrapolate voltage copies
+    EXTRAPOLATE_VOLTAGE = 0,
+
+    # Maximum number of fourier coefficients saved in file name
+    FOURIER_COEF_MAX = 3,
+
     # verbose logging for performance analysis
     VERBOSE = 0,
 
@@ -33,6 +45,9 @@ defaults = dict(
     # 0.5 corresponds to the old behaviour.
     LAZY_INVERSE_FACTOR = 0.25,
 
+    # disable plotting (for jobs on the cluster)
+    DISABLE_PLOTTING = 0,
+
     # use rtrg_cublas instead of rtrg_c
     USE_CUBLAS = 0,
 )