diff --git a/final_plots.py b/final_plots.py
index 837063880223a25af7c014fd30f9f6a88f0a5619..accc8ab3d8e67b96c8ba81e6c5d661ae001ff3aa 100644
--- a/final_plots.py
+++ b/final_plots.py
@@ -500,6 +500,238 @@ def export_interp(
)
+def export_interp_relative_o3a(
+ filename,
+ fixed_parameter,
+ fixed_value,
+ x_parameter,
+ y_parameter,
+ x_arr,
+ y_arr,
+ x_sort_parameter = None,
+ y_sort_parameter = None,
+ x_func = None,
+ y_func = None,
+ kxorder = 2,
+ kyorder = 2,
+ spline_s = 2e-4,
+ special_selection = None,
+ o2_scale_x = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
+ o2_scale_y = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
+ o2_scale_fixed = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
+ **parameters
+ ):
+ """
+ Export interpolated data. Parameters:
+
+ filename: string, file path
+ fixed_parameter: name of fixed parameter (e.g. "omega")
+ fixed_value: value of fixed parameter (e.g. 16.5372)
+ x_parameter: name of x axis parameter (e.g. "vdc")
+ y_parameter: name of y axis parameter (e.g. "vac")
+ x_sort_parameter: name of x axis parameter for sorting
+ y_sort_parameter: name of y axis parameter for sorting
+ x_func: function to obtain x value from data
+ y_func: function to obtain y value from data
+ kxorder: order of spline interpolation in x direction
+ kyorder: order of spline interpolation in y direction
+ special_selection: function for filtering input data. Arguments:
+ data, order, method, padding, **kwargs
+
+ Data are stored with units e=hbar=kB=Tkrg=1.
+ All results are exported as 2d arrays of the same shape.
+ The parameters vdc and vac are included as 2d arrays.
+ Results have the naming convention
+
+ {observable}_{method}_{order}{padding}{suffix}
+
+ observable:
+ gdc DC differential conductance, may have suffix "_ispline"
+ gac AC differential conductance, must have suffix "_ispline"
+ idc average current, may have suffix "_spline"
+ iac oscillating current, may have suffix "_spline"
+ phase AC phase
+ gamma Γ(E=0) (not an observable)
+
+ order:
+ o2 2nd order RG equations
+ o3 3rd order RG equations without Ga and with approximated integral
+ o3a 3rd order RG equations with approximated integral
+ o3p 3rd order RG equations with full integral
+
+ padding:
+ "" no Floquet matrix extrapolation
+ "_p" Floquet matrix extrapolation to mitigate truncation effects
+
+ method:
+ mu no unitary transformation, oscillating chemical potentials
+ J include oscillating voltage in coupling by unitary transformation
+
+ suffix
+ "" interpolation is done using griddata
+ "_spline" interpolation is done using spline
+ "_ispline" data is derived from spline interpolation for the
+ current (dc or ac)
+ """
+ assert isinstance(filename, str)
+ if x_func is None:
+ x_func = lambda data: data[x_parameter]
+ if y_func is None:
+ y_func = lambda data: data[y_parameter]
+ x_mesh, y_mesh = np.meshgrid(x_arr, y_arr)
+ mesh = (x_mesh, y_mesh)
+ dm = DataManager()
+ data = dm.list(**{fixed_parameter:fixed_value}, **parameters)
+ global_bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
+ | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
+ | DataManager.SOLVER_FLAGS["deleted"]
+ data = data.loc[data.solver_flags & global_bad_flags == 0]
+ data = data.sort_values([x_sort_parameter or x_parameter, y_sort_parameter or y_parameter])
+ data.vdc = np.round(data.vdc, 9)
+ data.vac = np.round(data.vac, 9)
+ required_flags = dict(
+ o3 = 0,
+ o3a = DataManager.SOLVER_FLAGS["include_Ga"],
+ o3p = DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
+ o3a_vb7 = DataManager.SOLVER_FLAGS["include_Ga"],
+ )
+ bad_flags = dict(
+ o3 = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
+ o3a = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
+ o3p = DataManager.SOLVER_FLAGS["second_order_rg_equations"],
+ o3a_vb7 = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
+ )
+ voltage_branches = dict(
+ o3 = 4,
+ o3a = 4,
+ o3p = 4,
+ o3a_vb7 = 7,
+ )
+ ref_data = data.loc[(data.solver_flags & (required_flags["o3a"] | bad_flags["o3a"]) == required_flags["o3a"]) & (data.method == "mu")]
+ merged = pd.merge(
+ data,
+ ref_data,
+ how="inner",
+ on=["vdc","vac","d","omega","energy_im","energy_re","solver_tol_rel","solver_tol_abs","xL","lazy_inverse_factor","resonant_dc_shift"],
+ suffixes=("", "_ref"))
+
+ def selection(order, method, padding=False, voltage_branches=4, **kwargs):
+ selected = (merged.solver_flags & required_flags[order] == required_flags[order]) \
+ & (merged.solver_flags & bad_flags[order] == 0) \
+ & (merged.method == method) \
+ & (merged.voltage_branches == voltage_branches)
+ if padding:
+ selected &= (merged.padding > 0) | (merged.nmax == 0)
+ else:
+ selected &= merged.padding == 0
+ for key, value in kwargs.items():
+ selected &= (merged[key] == value)
+ if special_selection is not None:
+ selected &= special_selection(merged, order, method, padding, **kwargs)
+ return merged.loc[selected]
+
+ def gen_data(order, method, padding=False, s=spline_s, extend_vdc=10, extend_vac=10, **kwargs):
+ suffix = "_p" if padding else ""
+ reduced_data = selection(order, method, padding, voltage_branches[order], **kwargs)
+ settings.logger.info(f"Starting with {order}{suffix} {method}, found {reduced_data.shape[0]} data points")
+ results = {}
+ if reduced_data.shape[0] < 100:
+ return results
+ if extend_vdc or extend_vac:
+ extended_data = [reduced_data]
+ if extend_vdc:
+ dc_mirrored = reduced_data[reduced_data.vdc<extend_vdc].copy()
+ dc_mirrored.vdc *= -1
+ dc_mirrored.dc_current *= -1
+ dc_mirrored.dc_current_ref *= -1
+ extended_data.append(dc_mirrored)
+ del dc_mirrored
+ if extend_vac:
+ ac_mirrored = reduced_data[reduced_data.vac<extend_vac].copy()
+ ac_mirrored.vac *= -1
+ ac_mirrored.ac_current_phase += np.pi
+ ac_mirrored.ac_current_phase_ref += np.pi
+ extended_data.append(ac_mirrored)
+ if extend_vdc:
+ acdc_mirrored = ac_mirrored[ac_mirrored.vdc<extend_vdc].copy()
+ acdc_mirrored.vdc *= -1
+ acdc_mirrored.dc_current *= -1
+ acdc_mirrored.dc_current_ref *= -1
+ extended_data.append(acdc_mirrored)
+ del acdc_mirrored
+ del ac_mirrored
+ reduced_data = pd.concat(extended_data)
+ xy_data = (x_func(reduced_data), y_func(reduced_data))
+ try:
+ results[f"gdc_{method}_{order}{suffix}"] = griddata(
+ xy_data,
+ reduced_data.dc_conductance - reduced_data.dc_conductance_ref,
+ mesh,
+ method="cubic")
+ except:
+ settings.logger.exception(f"in griddata interpolation of Gdc for {order}{suffix} {method}")
+ try:
+ results[f"idc_{method}_{order}{suffix}"] = griddata(
+ xy_data,
+ reduced_data.dc_current - reduced_data.dc_current_ref,
+ mesh,
+ method="cubic")
+ except:
+ settings.logger.exception(f"in griddata interpolation of Idc for {order}{suffix} {method}")
+ try:
+ idc_tck = bisplrep(*xy_data[::-1], reduced_data.dc_current - reduced_data.dc_current_ref, s=s, kx=kxorder, ky=kyorder)
+ results[f"idc_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, idc_tck)
+ results[f"gdc_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, idc_tck, dy=1)
+ except:
+ settings.logger.exception(f"in spline interpolation of Idc for {order}{suffix} {method}")
+ try:
+ results[f"iac_{method}_{order}{suffix}"] = griddata(
+ xy_data,
+ 2*(reduced_data.ac_current_abs - reduced_data.ac_current_abs_ref),
+ mesh,
+ method="cubic")
+ except:
+ settings.logger.exception(f"in griddata interpolation of Iac for {order}{suffix} {method}")
+ try:
+ iac_tck = bisplrep(*xy_data[::-1], 2*(reduced_data.ac_current_abs - reduced_data.ac_current_abs_ref), s=s, kx=kxorder, ky=kyorder)
+ results[f"iac_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, iac_tck)
+ results[f"gac_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, iac_tck, dx=1)
+ except:
+ settings.logger.exception(f"in spline interpolation of Iac for {order}{suffix} {method}")
+ try:
+ sel = reduced_data.vac > 0
+ results[f"phase_{method}_{order}{suffix}"] = griddata(
+ (xy_data[0][sel], xy_data[1][sel]),
+ reduced_data.ac_current_phase[sel] - reduced_data.ac_current_phase_ref[sel],
+ mesh,
+ method="cubic")
+ except:
+ settings.logger.exception(f"in griddata interpolation of phase for {order}{suffix} {method}")
+ try:
+ results[f"gamma_{method}_{order}{suffix}"] = griddata(
+ xy_data,
+ reduced_data.gamma - reduced_data.gamma_ref,
+ mesh,
+ method="cubic")
+ except:
+ settings.logger.exception(f"in griddata interpolation of gamma for {order}{suffix} {method}")
+ return results
+
+ np.savez(
+ filename,
+ **{fixed_parameter:fixed_value, x_parameter:x_mesh, y_parameter:y_mesh},
+ **gen_data("o3", "mu"),
+ **gen_data("o3a_vb7", "mu"),
+ **gen_data("o3p", "mu", integral_method=-1),
+ **gen_data("o3", "J", False),
+ **gen_data("o3a", "J", False),
+ **gen_data("o3p", "J", False, integral_method=-1),
+ **gen_data("o3", "J", True),
+ **gen_data("o3a", "J", True),
+ **gen_data("o3p", "J", True, integral_method=-1),
+ )
+
+
def export_omega5_interp(
filename = "figdata/omega5_interp.npz",
omega = 16.5372,
@@ -521,7 +753,8 @@ def export_omega5_interp(
return (data["vdc"] < vdc_max_J + 1e-3) & (data["vac"] < vac_max_J + 1e-3)
return True
yarr = np.linspace(vac_min, vac_max, ac_res)
- yarr[-1] -= 1e-6
+ yarr[-1] -= 1e-10
+ yarr[ac_res//2] -= 1e-10
export_interp(filename,
fixed_parameter = "omega",
fixed_value = omega,
@@ -539,6 +772,44 @@ def export_omega5_interp(
xL = 0.5)
+def export_omega5_interp_deviation(
+ filename = "figdata/omega5_interp_deviation.npz",
+ omega = 16.5372,
+ vdc_min = 0,
+ vdc_max = 165.372,
+ vac_min = 0,
+ vac_max = 165.372,
+ dc_res = 301,
+ ac_res = 201,
+ vdc_max_J = 82.686,
+ vac_max_J = 82.686,
+ ):
+ """
+ Export interpolated data for fixed frequency omega.
+ """
+ def special_selection(data, order, method, padding, **kwargs):
+ if method == "J":
+ return (data["vdc"] < vdc_max_J + 1e-3) & (data["vac"] < vac_max_J + 1e-3)
+ return True
+ yarr = np.linspace(vac_min, vac_max, ac_res)
+ yarr[-1] -= 1e-10
+ yarr[ac_res//2] -= 1e-10
+ export_interp_relative_o3a(filename,
+ fixed_parameter = "omega",
+ fixed_value = omega,
+ x_parameter = "vdc",
+ y_parameter = "vac",
+ x_arr = np.linspace(vdc_min, vdc_max, dc_res),
+ y_arr = yarr,
+ kyorder = 3,
+ special_selection = special_selection,
+ spline_s = 2e-5,
+ d = 1e9,
+ solver_tol_rel = 1e-8,
+ solver_tol_abs = 1e-10,
+ xL = 0.5)
+
+
def export_vdc0_interp(
filename = "figdata/vdc0_interp.npz",
omega_min = 0.1,
@@ -604,11 +875,14 @@ def prepare_RGconvergence():
real_cmap = lambda x: plt.cm.seismic((x+1)/2)[...,:3]
imwrite(f"figdata/colorbar_seismic.png", np.array(0xffff*plt.cm.seismic(np.linspace(0, 1, 0x1000))[...,:3], dtype=np.uint16).reshape((1,0x1000,3)), format="PNG-FI")
+ imwrite(f"figdata/colorbar_seismic_vertical.png", np.array(0xffff*plt.cm.seismic(np.linspace(0, 1, 0x1000))[...,:3], dtype=np.uint16).reshape((0x1000,1,3)), format="PNG-FI")
def export_img(name, array):
imwrite(f"figdata/{name}.png", np.array(0xffff*real_cmap(array[::-1]), dtype=np.uint16), format="PNG-FI")
data = np.load("figdata/omega5_interp.npz")
data_vb7 = np.load("figdata/omega5_interp_vb7.npz")
+ data_rel = np.load("figdata/omega5_interp_deviation.npz")
+
# O2
omega = 16.5372
omega_o2 = 48.695054
@@ -627,13 +901,43 @@ def prepare_RGconvergence():
export_img(f"convergence_{method}_relative", grel/normalize)
print(f"Color scale for {method} relative: {normalize}")
- gdiff = data[f"gdc_mu_o3a"] - data_vb7[f"gdc_mu_o3a"]
+ for method, u, sign in (("o3p","mu",-1), ("o3a_vb7","mu",-1), ("o3a","J",1), ("o3a_p","J",1)):
+ gdiff = sign*data_rel[f"gdc_{u}_{method}"]
+ if u != "mu":
+ method += "_" + u
+ normalize = np.nanmax(np.abs(gdiff))
+ export_img(f"convergence_{method}_diff_fromdiff", gdiff/normalize)
+ print(f"Color scale for {method} difference (interpolate in the end): {np.pi*normalize}")
+ grel = gdiff/data[f"gdc_mu_o3a"]
+ normalize = np.nanmax(np.abs(grel))
+ export_img(f"convergence_{method}_relative_fromdiff", grel/normalize)
+ print(f"Color scale for {method} relative (interpolate in the end): {normalize}")
+
+ export_img(f"convergence_o3a_J_relative_fromdiff_pnorm", data_rel[f"gdc_J_o3a"]/normalize)
+
+ g1rel = data_rel["gdc_J_o3a"] / data["gdc_mu_o3a"]
+ g2rel = data_rel["gdc_J_o3a_p"] / data["gdc_mu_o3a"]
+ g2rel[:16,:8] = np.nan
+ normalize = max(np.nanmax(np.abs(g1rel)), np.nanmax(np.abs(g2rel)))
+ g1img = np.array(0xffff*real_cmap(g1rel[::-1]/normalize), dtype=np.uint16)
+ g2img = np.array(0xffff*real_cmap(g2rel[::-1]/normalize), dtype=np.uint16)
+ g2img[-16:,:8] = 0x8000
+ imwrite(f"figdata/convergence_o3a_J_relative_fromdiff_maskednorm.png", g1img, format="PNG-FI")
+ imwrite(f"figdata/convergence_o3a_J_p_relative_fromdiff_masked.png", g2img, format="PNG-FI")
+ print(f"Color scale for o3a J relative (interpolate in the end, masked): {normalize}")
+
+ grel = -data_rel["gdc_mu_o3a_vb7"]/data["gdc_mu_o3a"]
+ normalize = max(np.nanmax(np.abs(grel[16:])), np.nanmax(np.abs(grel[:16,8:])))
+ export_img(f"convergence_o3a_vb7_relative_fromdiff_masked", grel/normalize)
+ print(f"Color scale for vb7 masked relative (interpolate in the end): {normalize}")
+
+ gdiff = data["gdc_mu_o3a"] - data_vb7["gdc_mu_o3a"]
normalize = np.nanmax(np.abs(gdiff))
- export_img(f"convergence_o3a_vb7_diff", gdiff/normalize)
+ export_img("convergence_o3a_vb7_diff", gdiff/normalize)
print(f"Color scale for o3a_vb7 difference: {np.pi*normalize}")
- grel = gdiff/data_vb7[f"gdc_mu_o3a"]
+ grel = gdiff/data_vb7["gdc_mu_o3a"]
normalize = np.nanmax(np.abs(grel))
- export_img(f"convergence_o3a_vb7_relative", grel/normalize)
+ export_img("convergence_o3a_vb7_relative", grel/normalize)
print(f"Color scale for o3a_vb7 relative: {normalize}")
for method, u in (("o2","mu"), ("o3","mu"), ("o3a","mu"), ("o3p","mu"), ("o3a_ispline","mu")):
@@ -1155,6 +1459,7 @@ def export_pulse_charge_pgfplots(omega=1.5):
_h5files = set()
parameters = dict(
d = 1e9,
+ omega = omega,
method = "J",
voltage_branches = 0,
include_Ga = True,
@@ -1168,9 +1473,11 @@ def export_pulse_charge_pgfplots(omega=1.5):
| DataManager.SOLVER_FLAGS["deleted"],
pulse_type = "gauss_symmetric",
)
+ parameters2 = parameters.copy()
+ parameters2.update(omega=2, solver_tol_rel=2e-8, solver_tol_abs=2e-10)
for (row, kondo) in itertools.chain(
- load_all_pulses_full(dm, omega=omega, **parameters),
- load_all_pulses_full(dm, omega=1.9, **parameters),
+ load_all_pulses_full(dm, **parameters),
+ load_all_pulses_full(dm, **parameters2),
):
ref_time = min(4*row.pulse_duration*row.omega/(2*np.pi), 0.15)
nmax = kondo.nmax
@@ -1510,6 +1817,99 @@ def export_floquet_matrices_pgfplots():
export_img("floquet_matrix_g201_J_diff", g201_matrices[("J",False,False)] - g201_ref_j)
+def export_Ga_max_pgfplots():
+ """
+ Show relevance of Ga: Plot maximum matrix element of Ga as function
+ of Vdc and Vac for fixed Ω.
+ """
+ dm = DataManager()
+ omega = 16.5372
+ parameters = dict(
+ omega = omega,
+ d = 1e9,
+ include_Ga = True,
+ integral_method = -1,
+ solver_tol_rel = 1e-8,
+ solver_tol_abs = 1e-10,
+ xL = 0.5,
+ bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
+ | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
+ | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
+ | DataManager.SOLVER_FLAGS["deleted"],
+ good_flags = DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
+ | DataManager.SOLVER_FLAGS["include_Ga"]
+ )
+ vdc = []
+ vac = []
+ ga_max = []
+ g2_max = []
+ g200c = []
+ ga00_max_c = []
+ nmax_arr = []
+ argmax = []
+ for kondo in dm.load_all(**parameters):
+ try:
+ try:
+ nmax = kondo.nmax
+ ga = kondo.ga
+ g2 = kondo.g2
+ except AttributeError:
+ continue
+ vdc.append(kondo.vdc)
+ vac.append(kondo.vac)
+ nmax_arr.append(nmax)
+ if ga.ndim == 4:
+ idx = np.abs(ga).argmax()
+ ga_max.append(np.abs(ga.flat[idx]))
+ g2_max.append(np.abs(g2).max())
+ argmax.append(idx)
+ ga00_max_c.append(np.abs(ga[0,0,nmax-4:nmax+5,nmax-4:nmax+5]).max())
+ g200c.append(np.abs(g2[0,0,nmax,nmax]))
+ elif ga.ndim == 5:
+ vb = kondo.voltage_branches
+ idx = np.abs(ga[:,:,vb]).argmax()
+ ga_max.append(np.abs(ga[:,:,vb].flat[idx]))
+ g2_max.append(np.abs(g2[:,:,vb]).max())
+ argmax.append(idx)
+ ga00_max_c.append(np.abs(ga[0,0,vb,nmax-4:nmax+5,nmax-4:nmax+5]).max())
+ g200c.append(np.abs(g2[0,0,vb,nmax,nmax]))
+ else:
+ raise ValueError("Invalid shape: %s"%ga.shape)
+ except:
+ settings.logger.exception("Error while reading data:")
+
+ vdc = np.array(vdc)
+ vac = np.array(vac)
+ nmax = np.array(nmax_arr)
+ argmax = np.array(argmax)
+ ga_max = np.array(ga_max)
+ g2_max = np.array(g2_max)
+ g200c = np.array(g200c)
+
+ vdc_max = 165.372
+ vac_max = 165.372
+ vdc_num = 301
+ vac_num = 201
+ vdc_arr = np.linspace(0, vdc_max, vdc_num)
+ vac_arr = np.linspace(0, vac_max, vac_num)
+ vdc_mesh, vac_mesh = np.meshgrid(vdc_arr, vac_arr)
+ extent = (-0.5*vdc_max/(vdc_num-1), vdc_max*(1+0.5/(vdc_num-1)), -0.5*vac_max/(vac_num-1), vac_max*(1+0.5/(vac_num-1)))
+ grid_ga_max = griddata((vdc, vac), ga_max, (vdc_mesh, vac_mesh), method="cubic")
+ grid_g2_max = griddata((vdc, vac), g2_max, (vdc_mesh, vac_mesh), method="cubic")
+
+ def export_img(name, array):
+ imwrite(f"figdata/{name}.png", np.array(0xffff*plt.cm.viridis(array[::-1])[...,:3], dtype=np.uint16), format="PNG-FI")
+
+ print(f"Extent: xmin={extent[0]/omega:.6g}, xmax={extent[1]/omega:,.6g}, ymin={extent[2]/omega:.6g}, ymax={extent[3]/omega:,.6g}")
+ norm = np.nanmax(np.abs(grid_ga_max))
+ export_img("Ga_max_abs", grid_ga_max/norm)
+ print("Ga_max_abs norm:", norm)
+ grid_ga_max_rel = grid_ga_max/grid_g2_max**2
+ norm = np.nanmax(np.abs(grid_ga_max_rel))
+ export_img("Ga_max_rel", grid_ga_max_rel/norm)
+ print("Ga_max_rel norm:", norm)
+
+
if __name__ == "__main__":
diff --git a/gen_data.py b/gen_data.py
index 88ab362937326379f16810a08d8be840699cf0eb..3be70e36b197531e5061c0276982167c4357efc2 100644
--- a/gen_data.py
+++ b/gen_data.py
@@ -381,6 +381,9 @@ python -m frtrg_kondo.gen_data \\
settings.logger.error(f"Invalid value nmax={nmax}: must be ≥0 at Vdc={kondo_options['vdc']:.8g}, Vac={kondo_options['vac']:.8g}, Ω={kondo_options['omega']:.8g}")
continue
settings.logger.info(f"Starting with Vdc={kondo_options['vdc']:.8g}, Vac={kondo_options['vac']:.8g}, Ω={kondo_options['omega']:.8g}, nmax={kondo_options['nmax']}")
+ if kondo_options["padding"] < 0:
+ kondo_options["padding"] = int(kondo_options["nmax"]*0.667)
+ settings.logger.info(f"...using padding={kondo_options['padding']}")
settings.logger.debug(kondo_options)
settings.logger.debug(solver_options)
kondo = Kondo(**kondo_options)
@@ -429,6 +432,9 @@ def save_data_mp(queue, lock, solver_options, filename, include='all', find_pole
settings.logger.error(f"Invalid value nmax={kondo_options['nmax']}: must be ≥0 at Vdc={kondo_options['vdc']:.8g}, Vac={kondo_options['vac']:.8g}, Ω={kondo_options['omega']:.8g}")
continue
settings.logger.info(f"Starting with Vdc={kondo_options['vdc']:.8g}, Vac={kondo_options['vac']:.8g}, Ω={kondo_options['omega']:.8g}, nmax={kondo_options['nmax']}")
+ if kondo_options["padding"] < 0:
+ kondo_options["padding"] = int(kondo_options["nmax"]*0.667)
+ settings.logger.info(f"...using padding={kondo_options['padding']}")
settings.logger.debug(kondo_options)
settings.logger.debug(solver_options)
kondo = Kondo(**kondo_options)
diff --git a/plot.py b/plot.py
index e1a302ebdbb0cb1a087b7eee11cd43a61139365a..a9a3ba384c92854b2e1cffe5ea89787479d03936 100644
--- a/plot.py
+++ b/plot.py
@@ -1882,12 +1882,12 @@ def plot_Ga_max(
argmax = []
for kondo in dm.load_all(omega=omega, method=method, **parameters):
try:
- vdc.append(kondo.vdc)
- vac.append(kondo.vac)
nmax = kondo.nmax
- nmax_arr.append(nmax)
ga = kondo.ga
g2 = kondo.g2
+ vdc.append(kondo.vdc)
+ vac.append(kondo.vac)
+ nmax_arr.append(nmax)
if ga.ndim == 4:
idx = np.abs(ga).argmax()
ga_max.append(np.abs(ga.flat[idx]))
diff --git a/tikz/commutator_RG2.tex b/tikz/commutator_RG2.tex
index 5491144dff2c08b18f97bac2a30bcafea7c2e6df..d53708a3c4cd6499684980ab5f8e05963cd31f4d 100644
--- a/tikz/commutator_RG2.tex
+++ b/tikz/commutator_RG2.tex
@@ -32,7 +32,7 @@
ytick = {0, 0.01, 0.02},
yticklabels = {$0$, $0.01$, $0.02$},
scaled ticks = false,
- ylabel = {$\|[\chi^{-1},ZG^2_{LR}]\|\,/\,J$},
+ ylabel = {$\|[\chi^{-1},ZG^2_{LR}]\|_\infty\,/\,J$},
]
\addplot graphics[xmin=0, xmax=1, ymin=0, ymax=0.0288735933] {../figdata/colorbar_viridis.png};
\end{axis}
diff --git a/tikz/convergence-j.tex b/tikz/convergence-j.tex
new file mode 100644
index 0000000000000000000000000000000000000000..1dcb2eaee95749365cd96f42c7992d28b08f2abd
--- /dev/null
+++ b/tikz/convergence-j.tex
@@ -0,0 +1,62 @@
+\begin{tikzpicture}
+ \pgfplotsset{set layers = axis on top}
+ \begin{axis}[
+ name = J,
+ anchor = west,
+ width = 4cm,
+ height = 4cm,
+ scale only axis,
+ title = {no truncation mitigation},
+ xmin = 0, xmax = 5,
+ ymin = 0, ymax = 5,
+ xlabel = {$\vdc~(\Omega)$},
+ ylabel = {$\vac~(\Omega)$},
+ clip = true,
+ ]
+ \addplot graphics[xmin=-.03333333, xmax=10.03333333, ymin=-.05, ymax=10.05] {../figdata/convergence_o3a_J_relative_fromdiff_maskednorm.png};
+ \node[anchor=north west] at (axis description cs:0.05,0.95) {(a)};
+ \end{axis}
+ \begin{axis}[
+ name = Jp,
+ anchor = west,
+ at = (J.east),
+ xshift = 4mm,
+ width = 4cm,
+ height = 4cm,
+ title = {with matrix extrapolation},
+ scale only axis,
+ xmin = 0, xmax = 5,
+ ymin = 0, ymax = 5,
+ xlabel = {$\vdc~(\Omega)$},
+ yticklabel = \empty,
+ clip = true,
+ ]
+ \addplot graphics[xmin=-.03333333, xmax=10.03333333, ymin=-.05, ymax=10.05] {../figdata/convergence_o3a_J_p_relative_fromdiff_masked.png};
+ \node[anchor=north west] at (axis description cs:0.05,0.95) {(b)};
+ \end{axis}
+ \begin{axis}[
+ name = colorbar,
+ at = (Jp.east),
+ anchor = west,
+ ymin = -0.00043767481879499436,
+ ymax = 0.00043767481879499436,
+ xshift = 8mm,
+ enlargelimits = false,
+ axis on top,
+ width = .75em,
+ height = 4cm,
+ scale only axis,
+ tickpos = right,
+ tick align = outside,
+ xtick = \empty,
+ %ytick = {-0.0004, 0, 0.0004},
+ %yticklabels = {$-0.04\%$, $0$, $0.04\%$},
+ %ytick = {-4e-4, -2e-4, 0, 2e-4, 4e-4},
+ %yticklabels = {$-0.04\%$, $-0.02\%$, $0$, $0.02\%$, $0.04\%$},
+ scaled ticks = false,
+ ylabel = {relative deviation},
+ %ylabel shift = -5mm,
+ ]
+ \addplot graphics[xmin=0, xmax=1, ymin=-0.00043767481879499436, ymax=0.00043767481879499436] {../figdata/colorbar_seismic_vertical.png};
+ \end{axis}
+\end{tikzpicture}
diff --git a/tikz/convergence-o2.tex b/tikz/convergence-o2.tex
index 22d27023903325761bd39a279ccbcfe0b9ac87ba..0b96d0a376b7e2932675506a5d299d35a88e8fd6 100644
--- a/tikz/convergence-o2.tex
+++ b/tikz/convergence-o2.tex
@@ -18,6 +18,7 @@
clip = false,
]
\addplot graphics[xmin=-.03333333, xmax=10.03333333, ymin=-.05, ymax=10.05] {../figdata/convergence_o2_diff.png};
+ \node at (2,9.2) {(b)};
% vertical lines
\addplot[no marks, black] coordinates {(0,1.6) (10,1.6)};
\addplot[no marks, black] coordinates {(0,2.6) (10,2.6)};
@@ -71,6 +72,7 @@
fill = none,
},
]
+ \node at (2,0.79) {(a)};
\addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o2_ac.dat};
\addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3p_ac.dat};
\addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o2_ac.dat};
@@ -108,6 +110,7 @@
},
clip = false,
]
+ \node at (1.3,9.2) {(c)};
\addplot[no marks, red] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
\addplot[no marks, red] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
\addplot[no marks, red] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o2_dc.dat};
diff --git a/tikz/convergence-o3.tex b/tikz/convergence-o3.tex
index e7484ad9fa60ca140c36714fa71dcc807e9083cf..280378b29bc78ab4fd37d0b8cf04ed750e90ff8d 100644
--- a/tikz/convergence-o3.tex
+++ b/tikz/convergence-o3.tex
@@ -18,6 +18,7 @@
clip = false,
]
\addplot graphics[xmin=-.03333333, xmax=10.03333333, ymin=-.05, ymax=10.05] {../figdata/convergence_o3_diff.png};
+ \node at (2,9.2) {(b)};
% vertical lines
\addplot[no marks, black] coordinates {(0,1.6) (10,1.6)};
\addplot[no marks, black] coordinates {(0,2.6) (10,2.6)};
@@ -71,6 +72,7 @@
fill = none,
},
]
+ \node at (2,0.79) {(a)};
\addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3_ac.dat};
\addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3p_ac.dat};
\addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o3_ac.dat};
@@ -108,6 +110,7 @@
},
clip = false,
]
+ \node at (1.3,9.2) {(c)};
\addplot[no marks, red] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
\addplot[no marks, red] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
\addplot[no marks, red] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o3_dc.dat};
diff --git a/tikz/convergence-o3a.tex b/tikz/convergence-o3a.tex
index bc277851d81decdfa593a8db1f4723b8b8d8c7a5..789f5770bf1682578d4183ddbf36029d0316ebaf 100644
--- a/tikz/convergence-o3a.tex
+++ b/tikz/convergence-o3a.tex
@@ -17,7 +17,8 @@
},
clip = false,
]
- \addplot graphics[xmin=-.03333333, xmax=10.03333333, ymin=-.05, ymax=10.05] {../figdata/convergence_o3a_diff.png};
+ \addplot graphics[xmin=-.03333333, xmax=10.03333333, ymin=-.05, ymax=10.05] {../figdata/convergence_o3p_diff_fromdiff.png};
+ \node at (2,9.2) {(b)};
% vertical lines
\addplot[no marks, black] coordinates {(0,1.6) (10,1.6)};
\addplot[no marks, black] coordinates {(0,2.6) (10,2.6)};
@@ -71,6 +72,7 @@
fill = none,
},
]
+ \node at (2,0.79) {(a)};
\addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3a_ac.dat};
\addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3p_ac.dat};
\addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o3a_ac.dat};
@@ -108,6 +110,7 @@
},
clip = false,
]
+ \node at (1.3,9.2) {(c)};
\addplot[no marks, red] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
\addplot[no marks, red] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
\addplot[no marks, red] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o3a_dc.dat};
@@ -130,8 +133,8 @@
\begin{axis}[
name = colorbar,
at = (dc.north west),
- xmin = -0.003211587341,
- xmax = 0.003211587341,
+ xmin = -0.002963806506594897,
+ xmax = 0.002963806506594897,
xshift = 5mm,
yshift = 11mm,
anchor = south west,
@@ -148,6 +151,6 @@
xlabel = {deviation $(2e^2/h)$},
xlabel shift = -15.85mm,
]
- \addplot graphics[ymin=0, ymax=1, xmin=-0.003211587341, xmax=0.003211587341] {../figdata/colorbar_seismic.png};
+ \addplot graphics[ymin=0, ymax=1, xmin=-0.002963806506594897, xmax=0.002963806506594897] {../figdata/colorbar_seismic.png};
\end{axis}
\end{tikzpicture}
diff --git a/tikz/convergence-o3a_ispline.tex b/tikz/convergence-o3a_ispline.tex
index 39f4d2156cf404a071e63847830ed69f4bd175cc..fdf693a06c42ee11a09430aa2e870c456615ce2e 100644
--- a/tikz/convergence-o3a_ispline.tex
+++ b/tikz/convergence-o3a_ispline.tex
@@ -18,6 +18,7 @@
clip = false,
]
\addplot graphics[xmin=-.03333333, xmax=10.03333333, ymin=-.05, ymax=10.05] {../figdata/convergence_o3a_ispline_relative.png};
+ \node at (2,9.2) {(b)};
% vertical lines
\addplot[no marks, black] coordinates {(0,1.6) (10,1.6)};
\addplot[no marks, black] coordinates {(0,2.6) (10,2.6)};
@@ -71,6 +72,7 @@
fill = none,
},
]
+ \node at (2,0.79) {(a)};
\addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3a_ispline_ac.dat};
\addplot[no marks, densely dotted] table [x=vdc, y expr={\thisrow{gdc1.60}+0.48}] {../figdata/convergence_mu_o3a_ac.dat};
\addplot[no marks, red] table [x=vdc, y expr={\thisrow{gdc2.60}+0.36}] {../figdata/convergence_mu_o3a_ispline_ac.dat};
@@ -108,6 +110,7 @@
},
clip = false,
]
+ \node at (1.3,9.2) {(c)};
\addplot[no marks, red] table [x expr={\thisrow{gdc0.00}+0.6}, y=vac] {../figdata/convergence_mu_o3a_ispline_dc.dat};
\addplot[no marks, red] table [x expr={\thisrow{gdc0.30}+0.48}, y=vac] {../figdata/convergence_mu_o3a_ispline_dc.dat};
\addplot[no marks, red] table [x expr={\thisrow{gdc0.70}+0.36}, y=vac] {../figdata/convergence_mu_o3a_ispline_dc.dat};
diff --git a/tikz/ga.tex b/tikz/ga.tex
new file mode 100644
index 0000000000000000000000000000000000000000..ee7ed0e989e7b4dfb083d54f18c038e2a0d9711a
--- /dev/null
+++ b/tikz/ga.tex
@@ -0,0 +1,82 @@
+\begin{tikzpicture}
+ \pgfplotsset{set layers = axis on top}
+ \begin{axis}[
+ name = abs,
+ anchor = west,
+ width = 4cm,
+ height = 4cm,
+ scale only axis,
+ title = {absolute value $\displaystyle\|G^a\|_\infty$},
+ xmin = 0, xmax = 10,
+ ymin = 0, ymax = 10,
+ xlabel = {$\vdc~(\Omega)$},
+ ylabel = {$\vac~(\Omega)$},
+ clip = true,
+ ]
+ \addplot graphics[xmin=-.03333333, xmax=10.03333333, ymin=-.05, ymax=10.05] {../figdata/Ga_max_abs.png};
+ \node[white,anchor=north west] at (axis description cs:0.05,0.95) {(a)};
+ \end{axis}
+ \begin{axis}[
+ name = cb_abs,
+ at = (abs.east),
+ anchor = west,
+ ymin = 0,
+ ymax = 0.110224621160373,
+ xshift = 4mm,
+ enlargelimits = false,
+ axis on top,
+ width = .75em,
+ height = 4cm,
+ scale only axis,
+ tickpos = right,
+ tick align = outside,
+ xtick = \empty,
+ scaled ticks = false,
+ ytick = {0,0.05,0.1},
+ yticklabels = {$0$,$0.05$,$0.1$},
+ %ylabel = {relative deviation},
+ %ylabel shift = -5mm,
+ ]
+ \addplot graphics[xmin=0, xmax=1, ymin=0, ymax=0.110224621160373] {../figdata/colorbar_viridis.png};
+ \end{axis}
+ %
+ \begin{axis}[
+ name = rel,
+ anchor = west,
+ at = (cb_abs.east),
+ xshift = 28mm,
+ width = 4cm,
+ height = 4cm,
+ title = {relative to $J^2$: $\displaystyle\frac{\|G^a\|_\infty}{\|G^2\|^2_\infty}$},
+ scale only axis,
+ xmin = 0, xmax = 10,
+ ymin = 0, ymax = 10,
+ xlabel = {$\vdc~(\Omega)$},
+ ylabel = {$\vac~(\Omega)$},
+ clip = true,
+ ]
+ \addplot graphics[xmin=-.03333333, xmax=10.03333333, ymin=-.05, ymax=10.05] {../figdata/Ga_max_rel.png};
+ \node[white,anchor=north west] at (axis description cs:0.05,0.95) {(b)};
+ \end{axis}
+ \begin{axis}[
+ name = cb_rel,
+ at = (rel.east),
+ anchor = west,
+ ymin = 0,
+ ymax = 0.8214882051039829,
+ xshift = 32mm,
+ enlargelimits = false,
+ axis on top,
+ width = .75em,
+ height = 4cm,
+ scale only axis,
+ tickpos = right,
+ tick align = outside,
+ xtick = \empty,
+ scaled ticks = false,
+ %ylabel = {relative deviation},
+ %ylabel shift = -5mm,
+ ]
+ \addplot graphics[xmin=0, xmax=1, ymin=0, ymax=0.8214882051039829] {../figdata/colorbar_viridis.png};
+ \end{axis}
+\end{tikzpicture}
diff --git a/tikz/harmonic_current.tex b/tikz/harmonic_current.tex
index aa1a866d1d5aefeaf2e8de51dbd451e797efd0bf..d11b0a4cfa58317d07560a857aa6a98b88c3f954 100644
--- a/tikz/harmonic_current.tex
+++ b/tikz/harmonic_current.tex
@@ -14,6 +14,7 @@
scaled ticks = false,
scale only axis,
]
+ \node[anchor=north west] at (axis description cs:0.04,0.96) {(a)};
\addplot[magenta, no marks] table [x=t, y=i5] {../figdata/harmonic_current-omega1.dat};
\addplot[blue, no marks] table [x=t, y=i10] {../figdata/harmonic_current-omega1.dat};
\addplot[green, no marks] table [x=t, y=i20] {../figdata/harmonic_current-omega1.dat};
@@ -44,6 +45,7 @@
scaled ticks = false,
scale only axis,
]
+ \node[anchor=north west] at (axis description cs:0.05,0.96) {(b)};
\addplot[magenta] table [x=p5, y=i5] {../figdata/harmonic_current-omega1.dat};
\addplot[blue] table [x=p10, y=i10] {../figdata/harmonic_current-omega1.dat};
\addplot[green] table [x=p20, y=i20] {../figdata/harmonic_current-omega1.dat};
@@ -70,6 +72,7 @@
scaled ticks = false,
scale only axis,
]
+ \node[anchor=north west] at (axis description cs:0.05,0.95) {(c)};
\addplot[blue] table [x=t, y=i16] {../figdata/rlm_current.dat};
\node[black, anchor=south] at (axis description cs:0.5,0.05) {\parbox{2cm}{\centering resonant\\level\\model}};
\end{axis}