diff --git a/data_management.py b/data_management.py
index e28ad52822cb4b64c510667aafa262886e66750d..a8841d01b1e89e87d6c85446a63fe3a86812cb70 100644
--- a/data_management.py
+++ b/data_management.py
@@ -145,6 +145,11 @@ class KondoExport:
                 solver_flags |= DataManager.SOLVER_FLAGS["simplified_initial_conditions"]
         except AttributeError:
             pass
+        try:
+            if self.kondo.improved_initial_conditions:
+                solver_flags |= DataManager.SOLVER_FLAGS["improved_initial_conditions"]
+        except AttributeError:
+            pass
         try:
             if self.kondo.include_Ga:
                 solver_flags |= DataManager.SOLVER_FLAGS["include_Ga"]
@@ -655,6 +660,7 @@ class DataManager:
             second_order_rg_equations = 0x0400,
             solve_integral_exactly = 0x0800,
             include_Ga = 0x1000,
+            improved_initial_conditions = 0x2000,
             )
 
     def __init__(self):
diff --git a/gen_data.py b/gen_data.py
index b9abb0c9080d9362ca0ba8ac8d9e311f1f920724..be69e3e9358b61d16ea5b4d7e8e4e9ab87626372 100644
--- a/gen_data.py
+++ b/gen_data.py
@@ -241,8 +241,11 @@ python -m frtrg_kondo.gen_data \\
     method_group.add_argument("--method", type=str, required=True, choices=('J', 'mu'),
             help = "J: include all time dependence in coupling by unitary transformation.\nmu: describe time dependence by Floquet matrix for chemical potentials.")
     method_group.add_argument("--simplified_initial_conditions", metavar="bool",
-            type=bool, default=False,
+            type=int, default=0,
             help = "Set initial condition for gammaL to 0")
+    method_group.add_argument("--improved_initial_conditions", metavar="bool",
+            type=int, default=0,
+            help = "Set initial condition for dgammaL/dE to linear response estimate")
     method_group.add_argument("--include_Ga", metavar="int", type=int,
             default=1, choices=(0,1),
             help = "include vertex parameter Ga in RG equations")
diff --git a/kondo.py b/kondo.py
index 8d785b6e0639797e049cc07a43b3cfd90aa8b879..ee5eaef3e9c4474835629830e36b1a14debbda6b 100644
--- a/kondo.py
+++ b/kondo.py
@@ -462,6 +462,7 @@ class Kondo:
             xL : 'asymmetry factor' = 0.5,
             compact = 0,
             simplified_initial_conditions = False,
+            improved_initial_conditions = False,
             truncation_order = 3,
             **rg_properties):
         """
@@ -510,6 +511,8 @@ class Kondo:
         simplified_initial_conditions : use simplified initial
                 conditions for the current kernel which lead to the same
                 result in the limit of large D.
+        improved_initial_conditions : use nonzero initial condition for
+                derivative of current kernel Γ^L
         truncation_order : Truncation order of RG equations, must be 2 or 3.
         """
         self.global_properties = GlobalRGproperties(
@@ -528,6 +531,8 @@ class Kondo:
         self.global_settings = settings.export()
         self.unitary_transformation = unitary_transformation
         self.simplified_initial_conditions = simplified_initial_conditions
+        self.improved_initial_conditions = improved_initial_conditions
+        assert not (simplified_initial_conditions and improved_initial_conditions)
         self.include_Ga = include_Ga
         self.solve_integral_exactly = solve_integral_exactly
         if truncation_order == 2:
@@ -737,9 +742,18 @@ class Kondo:
                 np.zeros((2*self.nmax+1,2*self.nmax+1), dtype=np.complex128),
                 symmetry = -symmetry
                 )
+        mu_matrix = self.mu.reduced(shift=1)
+        if self.improved_initial_conditions:
+            gamma = gamma0[self.voltage_branches].diagonal()
+            j = j0[self.voltage_branches].diagonal()
+            z = z0[self.voltage_branches].diagonal()
+            ddGammaL = 12j*np.pi*self.xL*(1-self.xL) * j**3/(self.d+gamma)
+            if self.truncation_order >= 3:
+                ddGammaL *= z**2*(1-z*j)
+            self.yL = RGfunction(self.global_properties, mu_matrix.values * ddGammaL.reshape((1,-1)), symmetry=-symmetry)
 
         ### Full current, also includes AC current
-        self.gammaL = self.mu.reduced(shift=1) @ self.deltaGammaL
+        self.gammaL = mu_matrix @ self.deltaGammaL
         if self.simplified_initial_conditions:
             self.gammaL *= 0
 
@@ -791,6 +805,11 @@ class Kondo:
         zj0 = z0 * j0 if self.truncation_order >= 3 else j0
         zjvalues = zvalues * jvalues if self.truncation_order >= 3 else jvalues
 
+        gamma0_red = gamma0[self.voltage_branches] if self.voltage_branches else gamma0
+        z0_red = z0[self.voltage_branches] if self.voltage_branches else z0
+        j0_red = j0[self.voltage_branches] if self.voltage_branches else j0
+        zj0_red = zj0[self.voltage_branches] if self.voltage_branches else zj0
+
         RGclass = SymRGfunction if self.compact else RGfunction
 
         # Create Γ and Z with the just calculated initial values.
@@ -912,16 +931,10 @@ class Kondo:
         self.current[0,0].symmetry = -symmetry
         self.current[1,1] = self.current[0,0].copy()
         if self.vac or self.resonant_dc_shift or self.fourier_coef is not None:
-            if self.voltage_branches:
-                if self.truncation_order >= 3:
-                    i0 = j0[self.voltage_branches] * (1 - zj0[self.voltage_branches])
-                else:
-                    i0 = j0[self.voltage_branches]
+            if self.truncation_order >= 3:
+                i0 = j0_red * (1 - zj0_red)
             else:
-                if self.truncation_order >= 3:
-                    i0 = j0 * (1 - zj0)
-                else:
-                    i0 = j0
+                i0 = j0_red
             i_LR = np.einsum(
                     'ij,...j->...ij',
                     init_matrix,
@@ -968,17 +981,10 @@ class Kondo:
                 )
         if self.resonant_dc_shift:
             assert self.compact == 0
-            if self.voltage_branches:
-                self.deltaGammaL.values[diag_idx] = 3*np.pi*self.xL*(1-self.xL) \
-                        * zj0[self.voltage_branches,self.resonant_dc_shift:-self.resonant_dc_shift]**2
-            else:
-                self.deltaGammaL.values[diag_idx] = 3*np.pi*self.xL*(1-self.xL) \
-                        * zj0[self.resonant_dc_shift:-self.resonant_dc_shift]**2
+            self.deltaGammaL.values[diag_idx] = 3*np.pi*self.xL*(1-self.xL) \
+                    * zj0_red[self.resonant_dc_shift:-self.resonant_dc_shift]**2
         else:
-            if self.voltage_branches:
-                diag_values = 3*np.pi*self.xL*(1-self.xL) * zj0[self.voltage_branches]**2
-            else:
-                diag_values = 3*np.pi*self.xL*(1-self.xL) * zj0**2
+            diag_values = 3*np.pi*self.xL*(1-self.xL) * zj0_red**2
             if self.compact:
                 assert diag_values.dtype == np.complex128
                 self.deltaGammaL.submatrix00 = np.diag(diag_values[0::2])
@@ -1018,19 +1024,39 @@ class Kondo:
                     mu[np.arange(i, 2*self.nmax+1), np.arange(2*self.nmax+1-i)] = f.conjugate()
             mu = RGclass(self.global_properties, mu, symmetry=1)
             self.gammaL = mu @ self.deltaGammaL.reduced()
+            if self.improved_initial_conditions:
+                ddGammaL = 12j*np.pi*self.xL*(1-self.xL) * j0_red**3/(self.d+gamma0_red)
+                if self.truncation_order >= 3:
+                    ddGammaL *= z0_red**2*(1-z0_red*j0_red)
+                if self.compact:
+                    self.yL = SymRGfunction(self.global_properties, mu.values * ddGammaL.reshape((1,-1)), symmetry=-symmetry, diag=False, offdiag=True)
+                else:
+                    self.yL = RGfunction(self.global_properties, mu.values * ddGammaL.reshape((1,-1)), symmetry=-symmetry)
         else:
             self.gammaL = self.vdc * self.deltaGammaL.reduced()
-            if self.vac and self.fourier_coef is None:
-                if self.voltage_branches:
-                    gammaL_AC = 3*np.pi*self.xL*(1-self.xL) * self.vac/2 * zj0[self.voltage_branches]**2
+            if self.improved_initial_conditions and self.vdc:
+                if self.truncation_order >= 3:
+                    self.yL.values += np.diag(self.vdc * 12j*np.pi*self.xL*(1-self.xL) * j0_red**3*z0_red**2*(1-z0_red*j0_red)/(self.d+gamma0_red))
                 else:
-                    gammaL_AC = 3*np.pi*self.xL*(1-self.xL) * self.vac/2 * zj0**2
+                    self.yL.values += np.diag(self.vdc * 12j*np.pi*self.xL*(1-self.xL) * j0_red**3/(self.d+gamma0_red))
+            if self.vac and self.fourier_coef is None:
+                gammaL_AC = 3*np.pi*self.xL*(1-self.xL) * self.vac/2 * zj0_red**2
                 if self.resonant_dc_shift:
                     gammaL_AC = gammaL_AC[...,self.resonant_dc_shift:-self.resonant_dc_shift]
                 idx = (np.arange(1, 2*self.nmax+1), np.arange(2*self.nmax))
                 self.gammaL.values[idx] = gammaL_AC[...,1:]
                 idx = (np.arange(0, 2*self.nmax), np.arange(1, 2*self.nmax+1))
                 self.gammaL.values[idx] = gammaL_AC[...,:-1]
+                if self.improved_initial_conditions:
+                    yL_AC = 12j*np.pi*self.xL*(1-self.xL) * self.vac/2 * j0_red**3/(self.d+gamma0_red)
+                    if self.truncation_order >= 3:
+                        yL_AC *= z0_red**2*(1-z0_red*j0_red)
+                    if self.resonant_dc_shift:
+                        yL_AC = yL_AC[...,self.resonant_dc_shift:-self.resonant_dc_shift]
+                    idx = (np.arange(1, 2*self.nmax+1), np.arange(2*self.nmax))
+                    self.yL.values[idx] = yL_AC[...,1:]
+                    idx = (np.arange(0, 2*self.nmax), np.arange(1, 2*self.nmax+1))
+                    self.yL.values[idx] = yL_AC[...,:-1]
         if self.simplified_initial_conditions:
             self.gammaL *= 0
 
diff --git a/settings.py b/settings.py
index 81f8e9d9c4e5fab475e391f061a9d069bda2ecdb..f9d8017242a4668074eda7f8ac7322c8135ad872 100644
--- a/settings.py
+++ b/settings.py
@@ -86,7 +86,7 @@ class GlobalFlags:
         BASEPATH = os.path.abspath("data"),
         DB_CONNECTION_STRING = "sqlite:///" + os.path.join(os.path.abspath("data"), "frtrg.sqlite"),
         FILENAME = "frtrg-04.h5",
-        VERSION = (14, 15, -1, -1),
+        VERSION = (14, 16, -1, -1),
         MIN_VERSION = (14, 0),
         LOG_TIME = 10, # in s
         ENFORCE_SYMMETRIC = 0,