.. _sphx_glr_dimer_lattice_plot_re_T_dos.py: =========================================== Evolution of DOS as function of temperature =========================================== Using a real frequency solver in the IPT scheme the Density of states is tracked through the first orders transition. .. code-block:: python # Created Tue Jun 14 15:44:38 2016 # Author: Óscar Nájera from __future__ import division, absolute_import, print_function import numpy as np import scipy.signal as signal from scipy.integrate import trapz import matplotlib.pyplot as plt import dmft.common as gf import dmft.ipt_real as ipt from dmft.utils import optical_conductivity import slaveparticles.quantum.dos as dos plt.matplotlib.rcParams.update({'axes.labelsize': 22, 'xtick.labelsize': 14, 'ytick.labelsize': 14, 'axes.titlesize': 22}) def loop_beta(u_int, tp, betarange, seed='mott gap'): """Solves IPT dimer and return Im Sigma_AA, Re Simga_AB returns list len(betarange) x 2 Sigma arrays """ s = [] g = [] w = np.linspace(-6, 6, 2**13) dw = w[1] - w[0] gss = gf.semi_circle_hiltrans(w + 5e-3j - tp - 1) gsa = gf.semi_circle_hiltrans(w + 5e-3j + tp + 1) for beta in betarange: print('U: ', u_int, 'tp: ', tp, 'Beta', beta) nfp = dos.fermi_dist(w, beta) (gss, gsa), (ss, sa) = ipt.dimer_dmft(u_int, tp, nfp, w, dw, gss, gsa) g.append((gss, gsa)) s.append((ss, sa)) return np.array(g), np.array(s) def simulation(U, tp, betarange): try: storage = np.load('dimer_Tevol_U{}tp{}.npz'.format(U, tp)) betarange = storage['betarange'] gwi = storage['gwi'] swi = storage['swi'] except FileNotFoundError: gwi, swi = loop_beta(U, tp, betarange) np.savez('dimer_Tevol_U{}tp{}.npz'.format(U, tp), betarange=betarange, gwi=gwi, swi=swi) return gwi, swi, betarange def plot_spectralfunc(gwi, betarange, yshift=False): plt.figure() shift = 0 for (gss, gsa), beta in zip(gwi, betarange): gloc = .5 * (gss + gsa) if yshift: shift = 100 / beta plt.plot(w, shift - gloc.imag / np.pi) plt.xlabel(r'$\omega$') plt.ylabel(r'$A(\omega)$') plt.xlim([-3, 3]) Polarization Bubble and optical conductivity -------------------------------------------- .. code-block:: python def optical_cond(ss, sa, tp, w, beta): E = np.linspace(-1, 1, 61) dos = np.exp(-2 * E**2) / np.sqrt(np.pi / 2) de = E[1] - E[0] dosde = (dos * de).reshape(-1, 1) nf = gf.fermi_dist(w, beta) lat_Aa = (-1 / np.add.outer(-E, w + tp - sa)).imag / np.pi lat_As = (-1 / np.add.outer(-E, w - tp - ss)).imag / np.pi a = optical_conductivity(lat_Aa, lat_Aa, nf, w, dosde) b = optical_conductivity(lat_As, lat_As, nf, w, dosde) c = optical_conductivity(lat_Aa, lat_As, nf, w, dosde) d = optical_conductivity(lat_As, lat_Aa, nf, w, dosde) return np.sum((a, b, c, d), axis=0) s-SNIM Experiment ----------------- .. code-block:: python def permitivity(real_sig, betarange, unit_weight): dw = w[1] - w[0] epsilon = [] for res, beta in zip(real_sig, betarange): sig = signal.hilbert(res, len(res) * 4)[:len(res)] ep = 1 + 1j * sig / \ (w - dw / 2) * unit_weight * 0.0074 # final unit conversion in SI epsilon.append(ep) return epsilon def plot_epsilon(epsilon, w): plt.figure() for ep in epsilon: plt.plot(w, ep.real, '-', w, ep.imag, ':') plt.xlabel(r'$\omega$') plt.ylabel(r'$\varepsilon$') def aal_ef(eps_sample, t, e_tip=-1300 + 960j, a=20): beta = (eps_sample - 1) / (eps_sample + 1) alpha = 4 * np.pi * a ** 3 * (e_tip - 1) / (e_tip + 2) return alpha * (1 + beta) / (1 - alpha * beta / (16 * np.pi * (3 * a + 2 * a * np.cos(t))**3)) Spectral function Dispersion .. code-block:: python # eps_k = np.linspace(-1., 1., 61) # ss, sa = swi[18] # lat_gfs = 1 / np.add.outer(-eps_k, w - tp + 5e-2j - ss) # lat_gfa = 1 / np.add.outer(-eps_k, w + tp + 5e-2j - sa) # Aw = np.clip(-.5 * (lat_gfa + lat_gfs).imag / np.pi, 0, 8) # plot_band_dispersion(w, Aw, 'Local ', eps_k, 'intensity') # plot_band_dispersion(w, -lat_gfa.imag / np.pi, 'Local ', eps_k, 'intensity') # plt.ylim([-3, 3]) .. code-block:: python w = np.linspace(-6, 6, 2**13) t = np.linspace(-np.pi, np.pi, 600) cos2 = np.cos(2 * t) # Reference insulator data from http://dx.doi.org/10.1126/science.1150124 s_ins = np.loadtxt( '/home/oscar/Dropbox/org/phd/dimer_lattice/optical conductivity/ins_340.csv', comments='x', delimiter=',').T s_exp_weight = trapz(s_ins[1], x=s_ins[0] / 8056.54) U = 2.8 tp = 0.3 higb = np.array([100., 70., 60., 50.]) temp = np.linspace(0.023, 0.035, 30) .. code-block:: python sig_wei = 20 .. code-block:: python U = 2.5 tp = 0.3 higb = np.array([100., 80., 70.]) temp = np.linspace(0.017, 0.021, 30) sig_wei = 3 temp = np.linspace(0.04, 0.05, 30) temp = np.linspace(0.023, 0.035, 30) .. code-block:: python betarange = np.concatenate((higb, 1 / temp, 1 / temp[::-1], higb[::-1])) gwi, swi, betarange = simulation(U, tp, betarange) resi = [optical_cond(ss, sa, tp, w, beta) for (ss, sa), beta in zip(swi, betarange)] Plots ----- .. code-block:: python temp_set = [0, 1, 2, 3, 4, len(temp)] plot_spectralfunc(gwi, betarange, yshift=True) plt.yticks(100 / betarange[temp_set], np.around(1 / betarange, 3)[temp_set]) plt.title(r'$U={}$ $t_\perp={}$'.format(U, tp)) plt.savefig('dimer_sig_Tevo_U{}tp{}_iptre_yshift.pdf'.format(U, tp)) plot_spectralfunc(gwi, betarange, yshift=False) plt.title(r'$U={}$ $t_\perp={}$'.format(U, tp)) plt.savefig('dimer_dos_Tevo_U{}tp{}_iptre_ontop.pdf'.format(U, tp)) plt.xlim([-.6, .6]) plt.ylim([0, 0.11]) plt.savefig('dimer_dos_Tevo_U{}tp{}_iptre_gapzoom.pdf'.format(U, tp)) plt.xlim([.6, 1.8]) plt.ylim([0.11, 0.45]) plt.savefig('dimer_dos_Tevo_U{}tp{}_iptre_qp_zoom.pdf'.format(U, tp)) rdw = (w > 0) * (w < s_ins[0][-1] / 8056.54) sim_weight = trapz(resi[sig_wei][rdw], x=w[rdw]) unit_weight = s_exp_weight / sim_weight plt.figure() for res, beta in zip(resi, betarange): plt.plot(w, unit_weight * res) plt.plot(w, unit_weight * resi[sig_wei], lw=2, label=r'$T={:.3}$'.format(1 / betarange[sig_wei])) plt.plot(w, unit_weight * resi[len(temp) + 1], lw=2, label=r'$T={:.3}$'.format(1 / betarange[len(temp) + 1])) plt.xlabel(r'$\omega$') plt.ylabel(r'$\Re \sigma(\omega)[\Omega cm]^{-1}$') plt.axvline(s_ins[0][-1] / 8056.54) plt.ylim([0, 8000]) plt.xlim([0, 2]) plt.legend(loc=0) plt.title(r'$U={}$ $t_\perp={}$'.format(U, tp)) plt.savefig('dimer_resigma_Tevo_U{}tp{}_iptre.pdf'.format(U, tp)) epsilon = permitivity(resi, betarange, unit_weight) plot_epsilon(epsilon, w) plt.axvline(w[4174]) plt.axvline(w[4247]) plt.ylim([-40, 110]) #plt.ylim([-1, 10]) plt.xlim([0.01, .6]) plt.title(r'$U={}$ $t_\perp={}$'.format(U, tp)) plt.savefig('dimer_epsilon_Tevo_U{}tp{}_iptre.pdf'.format(U, tp)) .. rst-class:: sphx-glr-horizontal * .. image:: /dimer_lattice/images/sphx_glr_plot_re_T_dos_001.png :scale: 47 * .. image:: /dimer_lattice/images/sphx_glr_plot_re_T_dos_002.png :scale: 47 * .. image:: /dimer_lattice/images/sphx_glr_plot_re_T_dos_003.png :scale: 47 * .. image:: /dimer_lattice/images/sphx_glr_plot_re_T_dos_004.png :scale: 47 .. code-block:: python plt.figure() epr = np.array([ep[4174] for ep in epsilon]) epp = np.array([ep[4247] for ep in epsilon]) plt.plot(1 / betarange, epr.real, 'g', label=r'$\omega=0.115$') plt.plot(1 / betarange, epr.imag, "g--") plt.plot(1 / betarange, epp.real, 'r', label=r'$\omega=0.222$') plt.plot(1 / betarange, epp.imag, "r--") plt.xlabel('T') plt.ylabel(r'$\varepsilon_r$') plt.legend(loc=0) plt.title(r'$U={}$ $t_\perp={}$'.format(U, tp)) plt.savefig('dimer_epsilon_Tevo_U{}tp{}_iptre_fixfreq.pdf'.format(U, tp)) .. image:: /dimer_lattice/images/sphx_glr_plot_re_T_dos_005.png :align: center .. code-block:: python # This normalization because this is an intensity measure so just to have # readable numbers sam = np.array([abs(trapz(aal_ef(ep, t) * cos2, t)) for ep in epr]) / 1e4 sap = np.array([abs(trapz(aal_ef(ep, t) * cos2, t, -474.41 + 248.61j)) for ep in epp]) / 1e4 plt.figure() plt.plot(1 / betarange, sam, 'x-', label=r'$\omega=0.115$') plt.plot(1 / betarange, sap, 'r', label=r'$\omega=0.222$') plt.xlabel('T') plt.ylabel('2$^{nd}$ harmonic amplitude') plt.title(r'$U={}$ $t_\perp={}$'.format(U, tp)) plt.legend(loc=0) plt.savefig('dimer_s2_Tevo_U{}tp{}_iptre_fixfreq.pdf'.format(U, tp)) .. image:: /dimer_lattice/images/sphx_glr_plot_re_T_dos_006.png :align: center **Total running time of the script:** ( 0 minutes 57.246 seconds) .. container:: sphx-glr-footer .. container:: sphx-glr-download :download:`Download Python source code: plot_re_T_dos.py ` .. container:: sphx-glr-download :download:`Download Jupyter notebook: plot_re_T_dos.ipynb ` .. rst-class:: sphx-glr-signature `Generated by Sphinx-Gallery `_