From 12bb8e4ffad81064969c4a28ffb22d900ba9ad46 Mon Sep 17 00:00:00 2001
From: Valentin Bruch <valentin.bruch@rwth-aachen.de>
Date: Mon, 31 Oct 2022 12:39:23 +0100
Subject: [PATCH] fixed bug in initial conditions for gammaL; implemented
 compact with Ga

---
 compact_rtrg.py |   2 +-
 kondo.py        | 239 ++++++++++++++++++++++++++++++------------------
 2 files changed, 153 insertions(+), 88 deletions(-)

diff --git a/compact_rtrg.py b/compact_rtrg.py
index d79472c..64a91b2 100644
--- a/compact_rtrg.py
+++ b/compact_rtrg.py
@@ -42,7 +42,7 @@ class SymRGfunction(RGfunction):
             self.submatrix00 = np.identity(self.nmax+1, dtype=np.complex128)
             self.submatrix11 = np.identity(self.nmax, dtype=np.complex128)
         elif values is not None:
-            if diag is None or offdiag is None or settings.CHECK_SYMMETRIES:
+            if diag is None or offdiag is None:
                 self.values = values
             else:
                 self.setValues(values, diag, offdiag)
diff --git a/kondo.py b/kondo.py
index 95ed97f..f8a35c2 100644
--- a/kondo.py
+++ b/kondo.py
@@ -794,14 +794,27 @@ class Kondo:
         RGclass = SymRGfunction if self.compact else RGfunction
 
         # Create Γ and Z with the just calculated initial values.
-        self.gamma = RGclass(self.global_properties, gammavalues, symmetry=0)
-        self.z = RGclass(self.global_properties, zvalues, symmetry=symmetry)
+        if self.compact:
+            self.gamma = SymRGfunction(self.global_properties, gammavalues, symmetry=symmetry, diag=True, offdiag=False)
+            self.gamma.check_symmetry()
+            self.z = SymRGfunction(self.global_properties, zvalues, symmetry=symmetry, diag=True, offdiag=False)
+            self.z.check_symmetry()
+        else:
+            self.gamma = RGfunction(self.global_properties, gammavalues, symmetry=symmetry)
+            self.gamma.check_symmetry()
+            self.z = RGfunction(self.global_properties, zvalues, symmetry=symmetry)
+            self.z.check_symmetry()
 
         if self.include_Ga:
             # Create Ga:
             if self.global_properties.symmetric:
-                self.ga_scalar = RGclass(self.global_properties, np.zeros_like(jvalues), symmetry=-symmetry)
+                settings.logger.debug("Creating Ga scalar")
+                if self.compact:
+                    self.ga_scalar = SymRGfunction(self.global_properties, np.zeros_like(jvalues), symmetry=-symmetry, diag=True, offdiag=False)
+                else:
+                    self.ga_scalar = RGfunction(self.global_properties, np.zeros_like(jvalues), symmetry=-symmetry)
             else:
+                settings.logger.debug("Creating Ga matrix")
                 self.ga = ReservoirMatrix(self.global_properties, symmetry=symmetry)
                 # TODO: check whether this works for RGclass with self.compact
                 self.ga[0,0] = RGclass(self.global_properties, np.zeros_like(jvalues), symmetry=-symmetry)
@@ -928,13 +941,24 @@ class Kondo:
             self.current[1,0] = -current_entry
 
         ## Initial conditions for voltage-variation of Γ: δΓ
-        self.deltaGamma = RGclass(
-                self.global_properties,
-                None if self.compact else np.zeros(
-                    (3,2*self.nmax+1,2*self.nmax+1) if self.voltage_branches else self.shape(),
-                    dtype=np.complex128),
-                symmetry = symmetry
-                )
+        if self.compact:
+            self.deltaGamma = SymRGfunction(
+                    self.global_properties,
+                    np.zeros(
+                        (3,2*self.nmax+1,2*self.nmax+1) if self.voltage_branches else self.shape(),
+                        dtype=np.complex128),
+                    symmetry = symmetry,
+                    diag = False,
+                    offdiag = True,
+                    )
+        else:
+            self.deltaGamma = RGfunction(
+                    self.global_properties,
+                    np.zeros(
+                        (3,2*self.nmax+1,2*self.nmax+1) if self.voltage_branches else self.shape(),
+                        dtype=np.complex128),
+                    symmetry = symmetry
+                    )
 
         ## Initial conditions for voltage-variation of current-Γ: δΓ_L
         self.deltaGammaL = RGclass(
@@ -965,11 +989,20 @@ class Kondo:
 
 
         ### Derivative of full current
-        self.yL = RGclass(
-                self.global_properties,
-                None if self.compact else np.zeros((2*self.nmax+1,2*self.nmax+1), dtype=np.complex128),
-                symmetry=-symmetry
-                )
+        if self.compact:
+            self.yL = SymRGfunction(
+                    self.global_properties,
+                    np.zeros((2*self.nmax+1,2*self.nmax+1), dtype=np.complex128),
+                    symmetry=-symmetry,
+                    diag=False,
+                    offdiag=True,
+                    )
+        else:
+            self.yL = RGfunction(
+                    self.global_properties,
+                    np.zeros((2*self.nmax+1,2*self.nmax+1), dtype=np.complex128),
+                    symmetry=-symmetry
+                    )
 
         ### Full current, also includes AC current
         if self.resonant_dc_shift == 0:
@@ -983,7 +1016,7 @@ class Kondo:
                 for i, f in enumerate(self.fourier_coef, 1):
                     mu[np.arange(2*self.nmax+1-i), np.arange(i, 2*self.nmax+1)] = f
                     mu[np.arange(i, 2*self.nmax+1), np.arange(2*self.nmax+1-i)] = f.conjugate()
-            mu = RGfunction(self.global_properties, mu, symmetry=-1)
+            mu = RGclass(self.global_properties, mu, symmetry=1)
             self.gammaL = mu @ self.deltaGammaL.reduced()
         else:
             self.gammaL = self.vdc * self.deltaGammaL.reduced()
@@ -1001,16 +1034,6 @@ class Kondo:
         if self.simplified_initial_conditions:
             self.gammaL *= 0
 
-        if self.compact:
-            self.deltaGamma.submatrix01 = np.zeros((self.nmax+1, self.nmax), dtype=np.complex128)
-            self.deltaGamma.submatrix10 = np.zeros((self.nmax, self.nmax+1), dtype=np.complex128)
-            self.yL.submatrix01 = np.zeros((self.nmax+1, self.nmax), dtype=np.complex128)
-            self.yL.submatrix10 = np.zeros((self.nmax, self.nmax+1), dtype=np.complex128)
-            self.gammaL.submatrix00 = None
-            self.gammaL.submatrix11 = None
-            self.gammaL.submatrix01 = np.zeros((self.nmax+1, self.nmax), dtype=np.complex128)
-            self.gammaL.submatrix10 = np.zeros((self.nmax, self.nmax+1), dtype=np.complex128)
-
         self.global_properties.energy = 1j*self.d
 
     def initialize(self, **solveopts):
@@ -1697,7 +1720,9 @@ class Kondo:
             shape11 = (nmax, nmax)
             shape01 = (nmax+1, nmax)
             shape10 = (nmax, nmax+1)
-            if self.include_Ga:
+            if self.include_Ga and hasattr(self, "ga_scalar"):
+                assert flattened_values.size == 8*(o+i) + 10*m + 6*f
+            elif self.include_Ga:
                 assert flattened_values.size == 9*(o+i) + 10*m + 8*f
             else:
                 assert flattened_values.size == 7*(o+i) + 10*m + 6*f
@@ -1730,20 +1755,18 @@ class Kondo:
             self.g3[0,1].values = flattened_values[7*(o+i)+10*m+2*f:7*(o+i)+10*m+3*f].reshape(shape_f)
             self.g3[1,0].values = flattened_values[7*(o+i)+10*m+3*f:7*(o+i)+10*m+4*f].reshape(shape_f)
             self.current[0,1].values = flattened_values[7*(o+i)+10*m+4*f:7*(o+i)+10*m+5*f].reshape(shape_f)
-            self.current[1,0].values = flattened_values[7*(o+i)+10*m+5*f:].reshape(shape_f)
+            self.current[1,0].values = flattened_values[7*(o+i)+10*m+5*f:7*(o+i)+10*m+6*f].reshape(shape_f)
             if self.include_Ga:
-                self.ga[0,0].submatrix00 = flattened_values[7*o+7*i+10*m+6*f:8*o+7*i+10*m+6*f].reshape(shape00)
-                self.ga[1,1].submatrix00 = flattened_values[8*o+7*i+10*m+6*f:9*o+7*i+10*m+6*f].reshape(shape00)
-                self.ga[0,0].submatrix11 = flattened_values[9*o+7*i+10*m+6*f:9*o+8*i+10*m+6*f].reshape(shape11)
-                self.ga[1,1].submatrix11 = flattened_values[9*o+8*i+10*m+6*f:9*o+9*i+10*m+6*f].reshape(shape11)
-                self.ga[0,1].values = flattened_values[9*o+9*i+10*m+6*f:9*o+9*i+10*m+7*f].reshape(shape_f)
-                self.ga[1,0].values = flattened_values[9*o+9*i+10*m+7*f:9*o+9*i+10*m+8*f].reshape(shape_f)
-            #self.g2[0,1].setValues(flattened_values[7*(o+i)+10*m:7*(o+i)+10*m+f].reshape(shape_f), True, True)
-            #self.g2[1,0].setValues(flattened_values[7*(o+i)+10*m+f:7*(o+i)+10*m+2*f].reshape(shape_f), True, True)
-            #self.g3[0,1].setValues(flattened_values[7*(o+i)+10*m+2*f:7*(o+i)+10*m+3*f].reshape(shape_f), True, True)
-            #self.g3[1,0].setValues(flattened_values[7*(o+i)+10*m+3*f:7*(o+i)+10*m+4*f].reshape(shape_f), True, True)
-            #self.current[0,1].setValues(flattened_values[7*(o+i)+10*m+4*f:7*(o+i)+10*m+5*f].reshape(shape_f), True, True)
-            #self.current[1,0].setValues(flattened_values[7*(o+i)+10*m+5*f:].reshape(shape_f), True, True)
+                if flattened_values.size == 8*(o+i) + 10*m + 6*f:
+                    self.ga_scalar.submatrix00 = flattened_values[7*o+7*i+10*m+6*f:8*o+7*i+10*m+6*f].reshape(shape00)
+                    self.ga_scalar.submatrix11 = flattened_values[8*o+7*i+10*m+6*f:8*o+8*i+10*m+6*f].reshape(shape11)
+                else:
+                    self.ga[0,0].submatrix00 = flattened_values[7*o+7*i+10*m+6*f:8*o+7*i+10*m+6*f].reshape(shape00)
+                    self.ga[1,1].submatrix00 = flattened_values[8*o+7*i+10*m+6*f:9*o+7*i+10*m+6*f].reshape(shape00)
+                    self.ga[0,0].submatrix11 = flattened_values[9*o+7*i+10*m+6*f:9*o+8*i+10*m+6*f].reshape(shape11)
+                    self.ga[1,1].submatrix11 = flattened_values[9*o+8*i+10*m+6*f:9*o+9*i+10*m+6*f].reshape(shape11)
+                    self.ga[0,1].values = flattened_values[9*o+9*i+10*m+6*f:9*o+9*i+10*m+7*f].reshape(shape_f)
+                    self.ga[1,0].values = flattened_values[9*o+9*i+10*m+7*f:9*o+9*i+10*m+8*f].reshape(shape_f)
         elif self.compact == 2:
             nmax = self.nmax
             f = (2*nmax+1)**2
@@ -1751,8 +1774,9 @@ class Kondo:
             i = (nmax**2 + 1) // 2
             m = nmax*(nmax+1) // 2
             if self.include_Ga:
-                raise NotImplementedError
-            assert flattened_values.size == 5*(o+i) + 8*m + 3*f
+                assert flattened_values.size == 6*(o+i) + 8*m + 3*f
+            else:
+                assert flattened_values.size == 5*(o+i) + 8*m + 3*f
             unpack01 = lambda flat, sym: np.concatenate((flat, sym*flat[::-1].conjugate()), axis=None).reshape((nmax+1,nmax))
             unpack10 = lambda flat, sym: np.concatenate((flat, sym*flat[::-1].conjugate()), axis=None).reshape((nmax,nmax+1))
             if nmax % 2:
@@ -1782,15 +1806,16 @@ class Kondo:
             self.current[0,0].submatrix01 = unpack01(flattened_values[5*(o+i)+6*m:5*(o+i)+7*m], -1)
             self.current[0,0].submatrix10 = unpack10(flattened_values[5*(o+i)+7*m:5*(o+i)+8*m], -1)
             self.current[1,1] = self.current[0,0].copy()
-            #self.g2[0,1].setValues(flattened_values[5*(o+i)+8*m:5*(o+i)+8*m+f].reshape((2*nmax+1, 2*nmax+1)), True, True)
-            #self.g3[0,1].setValues(flattened_values[5*(o+i)+8*m+f:5*(o+i)+8*m+2*f].reshape((2*nmax+1, 2*nmax+1)), True, True)
-            #self.current[0,1].setValues(flattened_values[5*(o+i)+8*m+2*f:].reshape((2*nmax+1, 2*nmax+1)), True, True)
             self.g2[0,1].values = flattened_values[5*(o+i)+8*m:5*(o+i)+8*m+f].reshape((2*nmax+1, 2*nmax+1))
             self.g2[1,0] = self.g2[0,1].floquetConjugate()
             self.g3[0,1].values = flattened_values[5*(o+i)+8*m+f:5*(o+i)+8*m+2*f].reshape((2*nmax+1, 2*nmax+1))
             self.g3[1,0] = -self.g3[0,1].floquetConjugate()
-            self.current[0,1].values = flattened_values[5*(o+i)+8*m+2*f:].reshape((2*nmax+1, 2*nmax+1))
+            self.current[0,1].values = flattened_values[5*(o+i)+8*m+2*f:5*(o+i)+8*m+3*f].reshape((2*nmax+1, 2*nmax+1))
             self.current[1,0] = -self.current[0,1].floquetConjugate()
+            if self.include_Ga:
+                self.ga_scalar.submatrix00 = unpack00(flattened_values[5*(o+i)+8*m+3*f:6*o+5*i+8*m+3*f], -1)
+                self.ga_scalar.submatrix11 = unpack11(flattened_values[6*o+5*i+8*m+3*f:6*(o+i)+8*m+3*f], -1)
+
 
         if settings.CHECK_SYMMETRIES:
             self.check_symmetry()
@@ -1983,18 +2008,22 @@ class Kondo:
                         self.currentE[1,0].values,
                     ]
             if self.include_Ga:
-                all_data += [
-                        self.gaE[0,0].submatrix00,
-                        self.gaE[1,1].submatrix00,
-                        self.gaE[0,0].submatrix11,
-                        self.gaE[1,1].submatrix11,
-                        self.gaE[0,1].values,
-                        self.gaE[1,0].values,
-                        ]
+                try:
+                    all_data += [
+                            self.ga_scalarE.submatrix00,
+                            self.ga_scalarE.submatrix11,
+                            ]
+                except AttributeError:
+                    all_data += [
+                            self.gaE[0,0].submatrix00,
+                            self.gaE[1,1].submatrix00,
+                            self.gaE[0,0].submatrix11,
+                            self.gaE[1,1].submatrix11,
+                            self.gaE[0,1].values,
+                            self.gaE[1,0].values,
+                            ]
             return np.concatenate(all_data, axis=None)
         elif self.compact == 2:
-            if self.include_Ga:
-                raise NotImplementedError
             if settings.CHECK_SYMMETRIES:
                 assert self.gammaE.symmetry == -1
                 assert self.gammaE.submatrix01 is None
@@ -2034,7 +2063,7 @@ class Kondo:
                     self.currentE[0,0].submatrix01 = dummy
                     self.currentE[0,0].submatrix10 = dummy
             pack_flat = lambda x, sym: (x.reshape(-1)[:(x.size+1)//2] + sym*x.reshape(-1)[:x.size//2-1:-1].conjugate())/2
-            return np.concatenate((
+            all_data = [
                         pack_flat(self.gammaE.submatrix00, -1),
                         pack_flat(self.zE.submatrix00, -1),
                         pack_flat(self.deltaGammaLE.submatrix00, -1),
@@ -2056,7 +2085,17 @@ class Kondo:
                         self.g2E[0,1].values,
                         self.g3E[0,1].values,
                         self.currentE[0,1].values,
-                ), axis=None)
+                    ]
+            if self.include_Ga:
+                assert self.ga_scalarE.submatrix01 is None
+                assert self.ga_scalarE.submatrix10 is None
+                assert np.allclose(self.ga_scalarE.submatrix00.flat, self.ga_scalarE.submatrix00.flat[::-1].conjugate())
+                assert np.allclose(self.ga_scalarE.submatrix11.flat, self.ga_scalarE.submatrix11.flat[::-1].conjugate())
+                all_data += [
+                        pack_flat(self.ga_scalarE.submatrix00, 1),
+                        pack_flat(self.ga_scalarE.submatrix11, 1),
+                        ]
+            return np.concatenate(all_data, axis=None)
 
 
     def packFlattenedDerivativesMinimal(self):
@@ -2102,14 +2141,20 @@ class Kondo:
                         self.g2E[1,0].values,
                     ]
             if self.include_Ga:
-                all_data += [
-                        self.gaE[0,0].submatrix00,
-                        self.gaE[1,1].submatrix00,
-                        self.gaE[0,0].submatrix11,
-                        self.gaE[1,1].submatrix11,
-                        self.gaE[0,1].values,
-                        self.gaE[1,0].values,
-                        ]
+                try:
+                    all_data += [
+                            self.ga_scalarE.submatrix00,
+                            self.ga_scalarE.submatrix11,
+                            ]
+                except AttributeError:
+                    all_data += [
+                            self.gaE[0,0].submatrix00,
+                            self.gaE[1,1].submatrix00,
+                            self.gaE[0,0].submatrix11,
+                            self.gaE[1,1].submatrix11,
+                            self.gaE[0,1].values,
+                            self.gaE[1,0].values,
+                            ]
             return np.concatenate(all_data, axis=None)
         elif self.compact == 2:
             if self.include_Ga:
@@ -2143,7 +2188,7 @@ class Kondo:
         (1d) array.
 
         Order of flattened_values:
-        Γ, Z, δΓ, *G2, *G3, *IL, δΓL, ΓL, YL
+        Γ, Z, δΓ, *G2, *G3, *IL, δΓL, ΓL, YL, [Ga]
         """
         if self.compact == 0:
             all_data = [
@@ -2222,18 +2267,22 @@ class Kondo:
                         self.current[1,0].values,
                     ]
             if self.include_Ga:
-                all_data += [
-                        self.ga[0,0].submatrix00,
-                        self.ga[1,1].submatrix00,
-                        self.ga[0,0].submatrix11,
-                        self.ga[1,1].submatrix11,
-                        self.ga[0,1].values,
-                        self.ga[1,0].values,
-                        ]
+                try:
+                    all_data += [
+                            self.ga_scalar.submatrix00,
+                            self.ga_scalar.submatrix11,
+                            ]
+                except AttributeError:
+                    all_data += [
+                            self.ga[0,0].submatrix00,
+                            self.ga[1,1].submatrix00,
+                            self.ga[0,0].submatrix11,
+                            self.ga[1,1].submatrix11,
+                            self.ga[0,1].values,
+                            self.ga[1,0].values,
+                            ]
             return np.concatenate(all_data, axis=None)
         elif self.compact == 2:
-            if self.include_Ga:
-                raise NotImplementedError
             if settings.CHECK_SYMMETRIES:
                 assert self.gamma.symmetry == 1
                 assert self.gamma.submatrix01 is None
@@ -2263,7 +2312,7 @@ class Kondo:
                 assert self.current[0,0].submatrix00 is None
                 assert self.current[0,0].submatrix11 is None
             pack_flat = lambda x, sym: (x.reshape(-1)[:(x.size+1)//2] + sym*x.reshape(-1)[:x.size//2-1:-1].conjugate())/2
-            return np.concatenate((
+            all_data = [
                         pack_flat(self.gamma.submatrix00, 1),
                         pack_flat(self.z.submatrix00, 1),
                         pack_flat(self.deltaGammaL.submatrix00, 1),
@@ -2285,7 +2334,17 @@ class Kondo:
                         self.g2[0,1].values,
                         self.g3[0,1].values,
                         self.current[0,1].values,
-                ), axis=None)
+                    ]
+            if self.include_Ga:
+                if settings.CHECK_SYMMETRIES:
+                    assert self.ga_scalar.symmetry == -1
+                    assert self.ga_scalar.submatrix01 is None
+                    assert self.ga_scalar.submatrix10 is None
+                all_data += [
+                        pack_flat(self.ga_scalar.submatrix00, -1),
+                        pack_flat(self.ga_scalar.submatrix11, -1),
+                        ]
+            return np.concatenate(all_data, axis=None)
 
 
     def packFlattenedValuesMinimal(self):
@@ -2332,14 +2391,20 @@ class Kondo:
                         self.g2[1,0].values,
                     ]
             if self.include_Ga:
-                all_data += [
-                        self.ga[0,0].submatrix00,
-                        self.ga[1,1].submatrix00,
-                        self.ga[0,0].submatrix11,
-                        self.ga[1,1].submatrix11,
-                        self.ga[0,1].values,
-                        self.ga[1,0].values,
-                        ]
+                try:
+                    all_data += [
+                            self.ga_scalar.submatrix00,
+                            self.ga_scalar.submatrix11,
+                            ]
+                except AttributeError:
+                    all_data += [
+                            self.ga[0,0].submatrix00,
+                            self.ga[1,1].submatrix00,
+                            self.ga[0,0].submatrix11,
+                            self.ga[1,1].submatrix11,
+                            self.ga[0,1].values,
+                            self.ga[1,0].values,
+                            ]
             return np.concatenate(all_data, axis=None)
         elif self.compact == 2:
             if self.include_Ga:
-- 
GitLab