diff --git a/final_plots.py b/final_plots.py index 2f1fba2a0171344e08954fdb7e277af81d552850..37d0d3a8ff9c3e61652e23659d877ad4b16aee77 100644 --- a/final_plots.py +++ b/final_plots.py @@ -862,6 +862,14 @@ def export_omega5_interp( xL = 0.5) +def export_omega5_interp_high_res(): + export_omega5_interp( + filename = "figdata/omega5_interp_high_res.npz", + dc_res = 601, + ac_res = 401, + ) + + def export_omega5_interp_deviation( filename = "figdata/omega5_interp_deviation.npz", omega = 16.5372, @@ -1059,19 +1067,23 @@ def prepare_RGconvergence(): if u != "mu": method += "_" + u g = np.pi*data[f"gdc_{u}_{method}"] - with open(f"figdata/convergence_{u}_{method}_dc.dat", "w") as file: - file.write("vac " + " ".join("gdc%.2f"%vdc[i] for i in dc_indices) + "\n") - np.savetxt(file, np.array([vac, *(g[:,i] for i in dc_indices)]).T) - with open(f"figdata/convergence_{u}_{method}_ac.dat", "w") as file: - file.write("vdc " + " ".join("gdc%.2f"%vac[i] for i in ac_indices) + "\n") - np.savetxt(file, np.array([vdc, *(g[i] for i in ac_indices)]).T) + np.savetxt(f"figdata/convergence_{u}_{method}_dc.dat", + np.array([vac, *(g[:,i] for i in dc_indices)]).T, + header="vac " + " ".join("gdc%.2f"%vdc[i] for i in dc_indices), + comments="") + np.savetxt(f"figdata/convergence_{u}_{method}_ac.dat", + np.array([vdc, *(g[i] for i in ac_indices)]).T, + header="vdc " + " ".join("gdc%.2f"%vac[i] for i in ac_indices), + comments="") g = np.pi*data_vb7[f"gdc_mu_o3a"] - with open(f"figdata/convergence_mu_o3a_vb7_dc.dat", "w") as file: - file.write("vac " + " ".join("gdc%.2f"%vdc[i] for i in dc_indices) + "\n") - np.savetxt(file, np.array([vac, *(g[:,i] for i in dc_indices)]).T) - with open(f"figdata/convergence_mu_o3a_vb7_ac.dat", "w") as file: - file.write("vdc " + " ".join("gdc%.2f"%vac[i] for i in ac_indices) + "\n") - np.savetxt(file, np.array([vdc, *(g[i] for i in ac_indices)]).T) + np.savetxt(f"figdata/convergence_mu_o3a_vb7_dc.dat", + np.array([vac, *(g[:,i] for i in dc_indices)]).T, + header="vac " + " ".join("gdc%.2f"%vdc[i] for i in dc_indices), + comments="") + np.savetxt(f"figdata/convergence_mu_o3a_vb7_ac.dat", + np.array([vdc, *(g[i] for i in ac_indices)]).T, + header="vdc " + " ".join("gdc%.2f"%vac[i] for i in ac_indices), + comments="") print("DC voltage: " + " ".join("%.6g"%vdc[i] for i in dc_indices)) print("AC voltage: " + " ".join("%.6g"%vac[i] for i in ac_indices)) @@ -1112,6 +1124,18 @@ def adjust_argument_array(start, end, *funcs, initial_resolution=51, s1=0.01, ma ys[i] = np.append(ys[i], func(new_x))[sorting] return (x, *ys) +def fourier(t, coef): + return 2*sum((c*np.exp(-2j*np.pi*t*n)).real for n,c in enumerate(coef)) - coef.real[0] + +def mklmbd(f, *args): + """ + Helper function to bind arguments to lambda function: + mklmbd(f, *args)(x) = f(x, *args) + """ + g = lambda x: f(x, *g.args) + g.args = args + return g + def export_gapproximations_pgfplots( omega = 16.5372, @@ -1130,7 +1154,9 @@ def export_gapproximations_pgfplots( | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \ | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \ | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \ - | DataManager.SOLVER_FLAGS["deleted"] + | DataManager.SOLVER_FLAGS["improved_initial_conditions"] \ + | DataManager.SOLVER_FLAGS["deleted"], + good_flags = DataManager.SOLVER_FLAGS["include_Ga"], ) dm = DataManager() @@ -1254,7 +1280,7 @@ def export_gapproximations_pgfplots( vdc_arr, g_pat_arr, gdc_arr, g_adiabatic_arr = adjust_argument_array( data_vdc.vdc.iloc[0], data_vdc.vdc.iloc[-1], g_pat_func, gdc_func, g_adiabatic_func, - **common_adjust_kw) + s1=2e-4, initial_resolution=101, maxit=5, scale_coef=0.9) np.savetxt("figdata/vdc_omega5_vac10.dat", np.array([vdc_arr/omega, gdc_arr, @@ -1304,6 +1330,7 @@ def export_asymmetry_pgfplots( | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \ | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \ | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \ + | DataManager.SOLVER_FLAGS["improved_initial_conditions"] \ | DataManager.SOLVER_FLAGS["deleted"] ).sort_values("vdc") data.drop_duplicates(("vdc", "xL"), inplace=True) @@ -1311,14 +1338,14 @@ def export_asymmetry_pgfplots( for x in xL: sel_data = data[np.abs(data.xL - x) < 1e-6] gdc_funcs.append(interp1d(sel_data.vdc, np.pi/(4*x*(1-x))*sel_data.dc_conductance, kind="cubic")) - vdc_arr, *gdc_arrs = adjust_argument_array(0, vdc_max, *gdc_funcs, s1=1e-3, yscales=[0.1 for x in xL], initial_resolution=101, maxit=3) + vdc_arr, *gdc_arrs = adjust_argument_array(0, vdc_max, *gdc_funcs, s1=2e-3, yscales=[0.1 for x in xL], initial_resolution=201, maxit=3) np.savetxt("figdata/asymmetry.dat", np.array([vdc_arr/omega, *gdc_arrs ]).T, header="vdc " + " ".join(f"xL{x:.3g}" for x in xL), comments = "", - fmt = "%.6g") + fmt = "%.9g") def export_kogan04_pgfplots(): @@ -1399,7 +1426,7 @@ def export_kogan04_pgfplots(): def export_bruhat18_pgfplots( plot_s = 2e-3, - init_res = 101, + init_res = 111, ): dm = DataManager() omega1 = 4.14 @@ -1423,6 +1450,7 @@ def export_bruhat18_pgfplots( | DataManager.SOLVER_FLAGS["second_order_rg_equations"] \ | DataManager.SOLVER_FLAGS["solve_integral_exactly"] \ | DataManager.SOLVER_FLAGS["extrapolate_voltage"] \ + | DataManager.SOLVER_FLAGS["improved_initial_conditions"] \ | DataManager.SOLVER_FLAGS["deleted"], ) gdc_funcs1 = [] @@ -1454,20 +1482,20 @@ def export_bruhat18_pgfplots( g1_pat_funcs = [get_pat(omega1, omega1*vac) for vac in vac_omega1] g2_pat_funcs = [get_pat(omega2, omega2*vac) for vac in vac_omega2] - vdc1, *arrays1 = adjust_argument_array(0, 21.5, *gdc_funcs1, s1=plot_s, maxit=3, initial_resolution=init_res) - vdc1_pat, *arrays1_pat = adjust_argument_array(0, 21.5, *g1_pat_funcs, s1=1.2*plot_s, maxit=3, initial_resolution=init_res) + vdc1, *arrays1 = adjust_argument_array(0, 22, *gdc_funcs1, s1=plot_s, maxit=3, initial_resolution=init_res) + vdc1_pat, *arrays1_pat = adjust_argument_array(0, 22, *g1_pat_funcs, s1=1.2*plot_s, maxit=3, initial_resolution=init_res) vdc1 = np.append(-vdc1[:0:-1], vdc1) / omega1 vdc1_pat = np.append(-vdc1_pat[:0:-1], vdc1_pat) / omega1 g1 = np.array([np.append(a[:0:-1], a) for a in arrays1]) g1_pat = np.array([np.append(a[:0:-1], a) for a in arrays1_pat]) - vdc2, *arrays2 = adjust_argument_array(0, 21.5, *gdc_funcs2, s1=plot_s, maxit=4, initial_resolution=init_res) - vdc2_pat, *arrays2_pat = adjust_argument_array(0, 21.5, *g2_pat_funcs, s1=1.2*plot_s, maxit=3, initial_resolution=init_res) + vdc2, *arrays2 = adjust_argument_array(0, 22, *gdc_funcs2, s1=plot_s, maxit=4, initial_resolution=init_res) + vdc2_pat, *arrays2_pat = adjust_argument_array(0, 22, *g2_pat_funcs, s1=1.2*plot_s, maxit=3, initial_resolution=init_res) vdc2 = np.append(-vdc2[:0:-1], vdc2) / omega2 vdc2_pat = np.append(-vdc2_pat[:0:-1], vdc2_pat) / omega2 g2 = np.array([np.append(a[:0:-1], a) for a in arrays2]) g2_pat = np.array([np.append(a[:0:-1], a) for a in arrays2_pat]) - vdc1_asym, *g1_asym = adjust_argument_array(0, 21.5, *gdc_funcs1_asym, s1=plot_s, maxit=3, initial_resolution=init_res) + vdc1_asym, *g1_asym = adjust_argument_array(0, 22, *gdc_funcs1_asym, s1=plot_s, maxit=3, initial_resolution=init_res) vdc1_asym = np.append(-vdc1_asym[:0:-1], vdc1_asym) / omega1 g1_asym = np.array([np.append(a[:0:-1], a) for a in g1_asym]) @@ -1478,7 +1506,7 @@ def export_bruhat18_pgfplots( # get background g_exp2_bg = g_exp2_20 - gdc_funcs2[0](np.abs(vdc_omega_exp2)*omega2) sel = ~np.isnan(g_exp2_bg) - g_exp2_bg_spl = splrep(vdc_omega_exp2[sel], g_exp2_bg[sel], s=1.1e-2) + g_exp2_bg_spl = splrep(vdc_omega_exp2[sel], g_exp2_bg[sel], s=1.15e-2) bg2 = np.array([quad_vec((lambda t: splev(vdc2+vac*np.cos(t), g_exp2_bg_spl, ext=3)), 0, np.pi)[0]/np.pi for vac in vac_omega2]) bg2_pat = np.array([quad_vec((lambda t: splev(vdc2_pat+vac*np.cos(t), g_exp2_bg_spl, ext=3)), 0, np.pi)[0]/np.pi for vac in vac_omega2]) bg1 = np.array([quad_vec((lambda t: splev((vdc1+vac*np.cos(t))*omega1/omega2, g_exp2_bg_spl, ext=3)), 0, np.pi)[0]/np.pi for vac in vac_omega1]) @@ -1564,13 +1592,26 @@ def export_pulse_current_pgfplots(omega=1.5, pulse_duration=0.01): np.linspace(0.5/t_prefactor, 3.33/t_prefactor, 40, endpoint=False), np.linspace(3.33/t_prefactor, 0.49, 41))) vdc, fourier_coef = fourier_coef_gauss_symmetric(1000, omega, pulse_duration, None, 1) + i_lmbds = [mklmbd(fourier, coef) for coef in i_coef] + g_lmbds = [mklmbd(fourier, coef) for coef in g_coef] + t1, *arrs1 = adjust_argument_array(-0.08/t_prefactor, 0.08/t_prefactor, *i_lmbds, *g_lmbds, s1=0.01, maxit=3, yscales=4*[0.0125]+4*[0.1]) + t2, *arrs2 = adjust_argument_array(0.08/t_prefactor, 0.5 - 0.08/t_prefactor, *i_lmbds, *g_lmbds, s1=1e-3, maxit=4, yscales=4*[1.5]+4*[1]) + t = np.append(t1, t2[1:]) + arrs = [] + for a1, a2 in zip(arrs1, arrs2): + arrs.append(np.append(a1, a2[1:])) + i_arrs = arrs[:4] + g_arrs = arrs[4:] + #i_arrs = [-coef[0].real + 2*sum((np.exp(-2j*np.pi*t*n)*c).real for n,c in enumerate(coef)) for coef in i_coef] + #g_arrs = [-coef[0].real + 2*sum((np.exp(-2j*np.pi*t*n)*c).real for n,c in enumerate(coef)) for coef in g_coef] u = vdc + 2*sum((np.exp(-2j*np.pi*t*n)*c).real for n,c in enumerate(fourier_coef, 1)) u /= u.max() - i_arrs = [-coef[0].real + 2*sum((np.exp(-2j*np.pi*t*n)*c).real for n,c in enumerate(coef)) for coef in i_coef] - g_arrs = [-coef[0].real + 2*sum((np.exp(-2j*np.pi*t*n)*c).real for n,c in enumerate(coef)) for coef in g_coef] t *= t_prefactor array = np.array([t, u, *i_arrs, *g_arrs]).T + start = np.searchsorted(t, -0.05, "right") - 1 + mid = np.searchsorted(t, 0.08, "right") + end = np.searchsorted(t, 3.33, "right") + 1 np.savetxt("figdata/pulse_current_full.dat", array, header = "t u " + " ".join(f"i{i}" for i in range(1,5)) + " " + " ".join(f"g{i}" for i in range(1,5)), @@ -1578,13 +1619,13 @@ def export_pulse_current_pgfplots(omega=1.5, pulse_duration=0.01): comments = "") np.savetxt("figdata/pulse_current_zoom.dat", - array[10:141], + array[start:mid+1], header = "t u " + " ".join(f"i{i}" for i in range(1,5)) + " " + " ".join(f"g{i}" for i in range(1,5)), fmt = "%.6g", comments = "") np.savetxt("figdata/pulse_current_tail.dat", - array[140:223], + array[mid-6:end], header = "t u " + " ".join(f"i{i}" for i in range(1,5)) + " " + " ".join(f"g{i}" for i in range(1,5)), fmt = "%.6g", comments = "") @@ -1703,13 +1744,13 @@ def export_pulse_charge_pgfplots(omega=1.5): def export_harmonic_modes( - #omega = 3.4425351, - omega = 68.850702, + omega = 3.4425351, + #omega = 68.850702, #omega = 16.5372, n_phase_crossings = 10, ): - label = "-omega20" - #label = "-omega1" + #label = "-omega20" + label = "-omega1" vac_omega = np.array([5, 10, 20, 40, 80]) dm = DataManager() data = dm.list( @@ -1729,13 +1770,13 @@ def export_harmonic_modes( good_flags = DataManager.SOLVER_FLAGS["improved_initial_conditions"] ) data = data.loc[data.method != "mu"] - current = [] - t = np.linspace(-0.5, 0.5, 201) - current = [] - phase = [] + #t = np.linspace(-0.5, 0.5, 201) + currents = [] + #phase = [] valid_vac_omega = [] phase_crossing_times = [] phase_crossing_currents = [] + fourier_scaled = lambda t, coef: (2*sum((c*np.exp(2j*np.pi*t*n)).real for n,c in enumerate(coef)) - coef.real[0]) / TK_VOLTAGE for vac in vac_omega: sel = np.abs(data.vac-vac*omega) < 1e-6 if sel.sum() == 0: @@ -1749,15 +1790,14 @@ def export_harmonic_modes( kondo, = KondoImport.read_from_h5(filename, row.hash) nmax = row.nmax coef = kondo.gammaL[nmax::-1,nmax] - iarr = 2*sum((c*np.exp(2j*np.pi*t*n)).real for n,c in enumerate(coef)) - coef.real[0] - iarr /= TK_VOLTAGE - current.append(iarr) + #currents.append((2*sum((c*np.exp(2j*np.pi*t*n)).real for n,c in enumerate(coef)) - coef.real[0]) / TK_VOLTAGE) + currents.append(mklmbd(fourier_scaled, kondo.gammaL[nmax::-1,nmax])) pc_times = np.arcsin(2*np.pi/vac*np.arange(n_phase_crossings)-1)/(2*np.pi) pc_currents = 2*sum((c*np.exp(2j*np.pi*pc_times*n)).real for n,c in enumerate(coef)) - coef.real[0] pc_currents /= TK_VOLTAGE phase_crossing_times.append(pc_times) phase_crossing_currents.append(pc_currents) - phase.append(vac/(2*np.pi)*(1+np.sin(2*np.pi*t))) + #phase.append(vac/(2*np.pi)*(1+np.sin(2*np.pi*t))) valid_vac_omega.append(vac) coef /= coef[1] np.savetxt( @@ -1766,6 +1806,12 @@ def export_harmonic_modes( header = "n i", comments = "", fmt = "%.6g") + t, *current = adjust_argument_array( + -0.5, 0.5, + *currents, + s1=0.4, initial_resolution=101, maxit=3, + yscales=(0.2,0.2)) + phase = [vac/(2*np.pi)*(1+np.sin(2*np.pi*t)) for vac in vac_omega] np.savetxt( f"figdata/harmonic_current{label}.dat", np.array([t, *current, *phase]).T, diff --git a/gen_data.py b/gen_data.py index da9b9ad14ed20a1b4a793497b172c1e5bcba215b..734aa3ef2049ddb410965c448c59294185cba7dc 100644 --- a/gen_data.py +++ b/gen_data.py @@ -250,7 +250,7 @@ python -m frtrg_kondo.gen_data \\ default=1, choices=(0,1), help = "include vertex parameter Ga in RG equations") method_group.add_argument("--solve_integral_exactly", metavar="bool", - type=bool, default=False, + type=int, default=0, help = "Solve integral in RG equations exactly by diagonalizing Floquet matrices. Requires --integral_method") method_group.add_argument("--integral_method", metavar="int", type=int, default=-15, help = "Select solution/approximation of frequency integral") @@ -285,20 +285,20 @@ python -m frtrg_kondo.gen_data \\ numerics_group.add_argument("--lazy_inverse_factor", metavar='float', type=float, default = settings.LAZY_INVERSE_FACTOR, help = "Factor between 0 and 1 for truncation of extended matrix before inversion.\n0 gives most precise results, 1 means discarding padding completely in inversion.\nOverwrites value set by environment variable LAZY_INVERSE_FACTOR.") - numerics_group.add_argument("--extrapolate_voltage", metavar='bool', type=bool, + numerics_group.add_argument("--extrapolate_voltage", metavar='bool', type=int, default = settings.EXTRAPOLATE_VOLTAGE, help = "Extrapolate along voltage branches (quadratic extrapolation).\nOverwrites value set by environment variable EXTRAPOLATE_VOLTAGE.") - numerics_group.add_argument("--check_symmetries", metavar='bool', type=bool, + numerics_group.add_argument("--check_symmetries", metavar='bool', type=int, default = settings.CHECK_SYMMETRIES, help = "Check symmetries during RG flow.\nOverwrites value set by environment variable CHECK_SYMMETRIES.") symmetry_group = numerics_group.add_mutually_exclusive_group() - symmetry_group.add_argument("--ignore_symmetries", metavar='bool', type=bool, + symmetry_group.add_argument("--ignore_symmetries", metavar='bool', type=int, default = settings.IGNORE_SYMMETRIES, help = "Do not use any symmetries.\nOverwrites value set by environment variable IGNORE_SYMMETRIES.") - symmetry_group.add_argument("--enforce_symmetric", metavar='bool', type=bool, + symmetry_group.add_argument("--enforce_symmetric", metavar='bool', type=int, default = settings.IGNORE_SYMMETRIES, help = "Enforce using symmetries, throw errors if no symmetries can be used.\nOverwrites value set by environment variable ENFORCE_SYMMETRIC.") - numerics_group.add_argument("--use_reference_implementation", metavar='bool', type=bool, + numerics_group.add_argument("--use_reference_implementation", metavar='bool', type=int, default = settings.USE_REFERENCE_IMPLEMENTATION, help = "Use slower reference implementation of RG equations instead of optimized implementation.\nOverwrites value set by environment variable USE_REFERENCE_IMPLEMENTATION.") diff --git a/kondo.py b/kondo.py index ee5eaef3e9c4474835629830e36b1a14debbda6b..93b7f7c1fd68a3456f83388a78cbcf1e22cf101b 100644 --- a/kondo.py +++ b/kondo.py @@ -209,11 +209,12 @@ def solveTV0_scalar( dj = dtheta*(1 + j/(1 + theta))/2 return np.array([dgamma, dtheta, dj]) + t_eval = solveopts.pop("t_eval", None) if full_output else (d,) result = solve_ivp( ode_function_imaxis, (0, d), np.array([gamma0, theta0, j0]), - t_eval = None if full_output else (d,), + t_eval = t_eval, rtol = rtol, atol = atol, **solveopts) @@ -1463,7 +1464,7 @@ class Kondo: # γ # dI γ ⎛ γ⊤ ⊤⎞⊤ 1 γ ⎛ ⎞ - # ——— = I Π G² + ⎜I Π G² ⎟ + — I ⎜ Π G² Z + Z G² Π ⎟ G³ + # ——— = I Π G² + ⎜I Π G² ⎟ + — I ⎜ Π G² Z + Z G² Π ⎟ G² # dE ⎝ ⎠ 2 34 ⎝ ⎠ 43 self.currentE = il @ pi @ g2 + ((il @ pi) % g2) + .5 * einsum_34_x_43(il, bracket, g2) @@ -1475,7 +1476,7 @@ class Kondo: # dE ⎝ ⎠ ⎝ ⎠ self.g2E += g2 @ pi @ ga + ga @ pi @ g2 - ((g2 @ pi) % ga) - ((ga @ pi) % g2) - # dG² ⎛ ⊤ ⊤⎞⊤ + # dG³ ⎛ ⊤ ⊤⎞⊤ # ——— += Ga Π G³ - ⎜Ga Π G³ ⎟ # dE ⎝ ⎠ self.g3E += ga @ pi @ g3 - ((ga @ pi) % g3) @@ -1536,12 +1537,6 @@ class Kondo: # ———— -= — tr ⎜I ⎜ Π δΓ Z + Z δΓ Π ⎟ G³⎟ # dE 4 ⎝ ⎝ ⎠ ⎠ self.deltaGammaLE += (-0.75) * (il @ (pi_reduced @ deltaGamma @ z_reduced + z_reduced @ deltaGamma @ pi_reduced) @ g3).tr() - elif self.truncation_order == 2: - # γ - # dδΓ 3 ⎛ ⎞ γ - # ———— = — i ⎜ δ - δ ⎟ I Π G³ - # dE 2 ⎝ 1L 2L ⎠ 12 21 - self.deltaGammaLE = 1.5j * (il[0,1] @ pi @ g3[1,0] - il[1,0] @ pi @ g3[0,1]) # γ # dΓ γ diff --git a/plot.py b/plot.py index af28de0ef82a738273ce4320cdb9003c750b1d2d..5b4f9997fbff053222254571940f601db7772ba7 100644 --- a/plot.py +++ b/plot.py @@ -1897,7 +1897,7 @@ def save_contour_coordinates(contour, filename, x, y, z): def RGeq_comparison_contours2( dm = None, - filename = "figdata/omega5_interp.npz", + filename = "figdata/omega5_interp_high_res.npz", observable = "gdc", prefactor = np.pi, **trashparams):