#+PROPERTY: header-args :session /ssh:l:/home/hiro/.local/share/jupyter/runtime/kernel-db283c80-f40c-4ded-8b78-99c9efe3be3c.json :kernel python :pandoc t :async yes * Setup ** Jupyter #+begin_src jupyter-python %load_ext autoreload %autoreload 2 %load_ext jupyter_spaces #+end_src #+RESULTS: : The autoreload extension is already loaded. To reload it, use: : %reload_ext autoreload ** Matplotlib #+begin_src jupyter-python import matplotlib import matplotlib.pyplot as plt #matplotlib.use("TkCairo", force=True) %gui tk %matplotlib inline plt.style.use('ggplot') #+end_src #+RESULTS: ** Richard (old) HOPS #+begin_src jupyter-python import hierarchyLib import hierarchyData import numpy as np from stocproc.stocproc import StocProc_FFT import bcf from dataclasses import dataclass, field import scipy import scipy.misc import scipy.signal import pickle from scipy.special import gamma as gamma_func from scipy.optimize import curve_fit #+end_src #+RESULTS: ** Auxiliary Definitions #+begin_src jupyter-python σ1 = np.matrix([[0,1],[1,0]]) σ2 = np.matrix([[0,-1j],[1j,0]]) σ3 = np.matrix([[1,0],[0,-1]]) #+end_src #+RESULTS: * Model Setup Basic parameters. #+begin_src jupyter-python class params: T = 2 t_max = 15 t_steps = int(t_max * 1/.05) k_max = 10 N = 4000 seed = 100 dim = 2 H_s = σ3 + np.eye(dim) L = σ2 #1 / 2 * (σ1 - 1j * σ2) ψ_0 = np.array([0, 1]) s = 1 num_exp_t = 4 wc = 1 with open("good_fit_data_abs_brute_force", "rb") as f: good_fit_data_abs = pickle.load(f) alpha = 0.8 # _, g_tilde, w_tilde = good_fit_data_abs[(numExpFit, s)] # g_tilde = np.array(g_tilde) # w_tilde = np.array(w_tilde) # g = 1 / np.pi * gamma_func(s + 1) * wc ** (s + 1) * np.asarray(g_tilde) # w = wc * np.asarray(w_tilde) bcf_scale = np.pi / 2 * alpha * wc ** (1 - s) #+end_src #+RESULTS: ** BCF and Thermal BCF #+begin_src jupyter-python @dataclass class CauchyBCF: δ: float wc: float def I(self, ω): return np.sqrt(self.δ) / (self.δ + (ω - self.wc) ** 2 / self.δ) def __call__(self, τ): return np.sqrt(self.δ) * np.exp(-1j * self.wc * τ - np.abs(τ) * self.δ) def __bfkey__(self): return self.δ, self.wc α = bcf.OBCF(s=params.s, eta=1, gamma=params.wc) I = bcf.OSD(s=params.s, eta=1, gamma=params.wc) #+end_src #+RESULTS: *** Fit We now fit a sum of exponentials against the BCF. #+begin_src jupyter-python from lmfit import minimize, Parameters def α_apprx(τ, g, w): return np.sum( g[np.newaxis, :] * np.exp(-w[np.newaxis, :] * (τ[:, np.newaxis])), axis=1 ) def _fit(): def residual(fit_params, x, data): resid = 0 w = np.array([fit_params[f"w{i}"] for i in range(params.num_exp_t)]) + 1j * np.array( [fit_params[f"wi{i}"] for i in range(params.num_exp_t)] ) g = np.array([fit_params[f"g{i}"] for i in range(params.num_exp_t)]) + 1j * np.array( [fit_params[f"gi{i}"] for i in range(params.num_exp_t)] ) resid = data - α_apprx(x, g, w) return resid.view(float) fit_params = Parameters() for i in range(params.num_exp_t): fit_params.add(f"g{i}", value=.1) fit_params.add(f"gi{i}", value=.1) fit_params.add(f"w{i}", value=.1) fit_params.add(f"wi{i}", value=.1) ts = np.linspace(0, params.t_max, 1000) out = minimize(residual, fit_params, args=(ts, α(ts))) w = np.array([out.params[f"w{i}"] for i in range(params.num_exp_t)]) + 1j * np.array( [out.params[f"wi{i}"] for i in range(params.num_exp_t)] ) g = np.array([out.params[f"g{i}"] for i in range(params.num_exp_t)]) + 1j * np.array( [out.params[f"gi{i}"] for i in range(params.num_exp_t)] ) return w, g w, g = _fit() #+end_src #+RESULTS: *** Plot Let's look a the result. #+begin_src jupyter-python class bcfplt: t = np.linspace(0, params.t_max, 1000) ω = np.linspace(params.wc - 10, params.wc + 10, 1000) fig, axs = plt.subplots(2) axs[0].plot(t, np.real(α(t))) axs[0].plot(t, np.imag(α(t))) axs[0].plot(t, np.real(α_apprx(t, g, w))) axs[0].plot(t, np.imag(α_apprx(t, g, w))) axs[1].plot(ω, I(ω).real) axs[1].plot(ω, I(ω).imag) #+end_src #+RESULTS: [[file:./.ob-jupyter/9f05a1fbf06920c271f0667db664ce2972415437.png]] Seems ok for now. ** Hops setup #+begin_src jupyter-python HierachyParam = hierarchyData.HiP( k_max=params.k_max, # g_scale=None, # sample_method='random', seed=params.seed, nonlinear=True, normalized=False, # terminator=False, result_type=hierarchyData.RESULT_TYPE_ALL, # accum_only=None, # rand_skip=None ) #+end_src #+RESULTS: Integration. #+begin_src jupyter-python IntegrationParam = hierarchyData.IntP( t_max=params.t_max, t_steps=params.t_steps, # integrator_name='zvode', # atol=1e-8, # rtol=1e-8, # order=5, # nsteps=5000, # method='bdf', # t_steps_skip=1 ) #+end_src #+RESULTS: And now the system. #+begin_src jupyter-python SystemParam = hierarchyData.SysP( H_sys=params.H_s, L=params.L, psi0=params.ψ_0, # excited qubit g=np.array(g), w=np.array(w), H_dynamic=[], bcf_scale=params.bcf_scale, # some coupling strength (scaling of the fit parameters 'g_i') gw_hash=None, # this is used to load g,w from some database len_gw=len(g), ) #+end_src #+RESULTS: The quantum noise. #+begin_src jupyter-python Eta = StocProc_FFT( I, params.t_max, α, negative_frequencies=False, seed=params.seed, intgr_tol=1e-3, intpl_tol=1e-3, scale=params.bcf_scale, ) #+end_src #+RESULTS: #+begin_example stocproc.stocproc - INFO - non neg freq only stocproc.method_ft - INFO - get_dt_for_accurate_interpolation, please wait ... stocproc.method_ft - INFO - acc interp N 33 dt 9.38e-01 -> diff 2.83e-01 stocproc.method_ft - INFO - acc interp N 65 dt 4.69e-01 -> diff 8.53e-02 stocproc.method_ft - INFO - acc interp N 129 dt 2.34e-01 -> diff 1.76e-02 stocproc.method_ft - INFO - acc interp N 257 dt 1.17e-01 -> diff 3.92e-03 stocproc.method_ft - INFO - acc interp N 513 dt 5.86e-02 -> diff 9.52e-04 stocproc.method_ft - INFO - requires dt < 5.859e-02 stocproc.method_ft - INFO - get_N_a_b_for_accurate_fourier_integral, please wait ... stocproc.method_ft - INFO - J_w_min:1.00e-02 N 32 yields: interval [0.00e+00,6.47e+00] diff 9.83e-03 stocproc.method_ft - INFO - J_w_min:1.00e-03 N 32 yields: interval [0.00e+00,9.12e+00] diff 8.12e-03 stocproc.method_ft - INFO - J_w_min:1.00e-02 N 64 yields: interval [0.00e+00,6.47e+00] diff 1.11e-02 stocproc.method_ft - INFO - increasing N while shrinking the interval does lower the error -> try next level stocproc.method_ft - INFO - J_w_min:1.00e-04 N 32 yields: interval [0.00e+00,1.17e+01] diff 1.32e-02 stocproc.method_ft - INFO - J_w_min:1.00e-03 N 64 yields: interval [0.00e+00,9.12e+00] diff 1.22e-03 stocproc.method_ft - INFO - J_w_min:1.00e-02 N 128 yields: interval [0.00e+00,6.47e+00] diff 1.14e-02 stocproc.method_ft - INFO - increasing N while shrinking the interval does lower the error -> try next level stocproc.method_ft - INFO - J_w_min:1.00e-05 N 32 yields: interval [0.00e+00,1.42e+01] diff 1.94e-02 stocproc.method_ft - INFO - J_w_min:1.00e-04 N 64 yields: interval [0.00e+00,1.17e+01] diff 2.57e-03 stocproc.method_ft - INFO - J_w_min:1.00e-03 N 128 yields: interval [0.00e+00,9.12e+00] diff 8.98e-04 stocproc.method_ft - INFO - return, cause tol of 0.001 was reached stocproc.method_ft - INFO - requires dx < 7.123e-02 stocproc.stocproc - INFO - Fourier Integral Boundaries: [0.000e+00, 1.251e+02] stocproc.stocproc - INFO - Number of Nodes : 2048 stocproc.stocproc - INFO - yields dx : 6.107e-02 stocproc.stocproc - INFO - yields dt : 5.023e-02 stocproc.stocproc - INFO - yields t_max : 1.028e+02 #+end_example The sample trajectories are smooth. #+begin_src jupyter-python %%space plot ts = np.linspace(0, params.t_max, 1000) Eta.new_process() plt.plot(ts, Eta(ts).real) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/10385a90753d7683d5740b88622cc4274de9c86a.png]] :END: And now the thermal noise. #+begin_src jupyter-python EtaTherm = StocProc_FFT( spectral_density=bcf.OFTDens(s=params.s, eta=1, gamma=params.wc, beta=1 / params.T), t_max=params.t_max, alpha=bcf.OFTCorr(s=params.s, eta=1, gamma=params.wc, beta=1 / params.T), intgr_tol=1e-3, intpl_tol=1e-3, seed=params.seed, negative_frequencies=False, scale=params.bcf_scale, ) #+end_src #+RESULTS: #+begin_example stocproc.stocproc - INFO - non neg freq only stocproc.method_ft - INFO - get_dt_for_accurate_interpolation, please wait ... stocproc.method_ft - INFO - acc interp N 33 dt 9.38e-01 -> diff 6.53e-02 stocproc.method_ft - INFO - acc interp N 65 dt 4.69e-01 -> diff 1.47e-02 stocproc.method_ft - INFO - acc interp N 129 dt 2.34e-01 -> diff 3.20e-03 stocproc.method_ft - INFO - acc interp N 257 dt 1.17e-01 -> diff 7.67e-04 stocproc.method_ft - INFO - requires dt < 1.172e-01 stocproc.method_ft - INFO - get_N_a_b_for_accurate_fourier_integral, please wait ... stocproc.method_ft - INFO - J_w_min:1.00e-02 N 32 yields: interval [0.00e+00,4.10e+00] diff 2.00e-02 stocproc.method_ft - INFO - J_w_min:1.00e-03 N 32 yields: interval [0.00e+00,5.82e+00] diff 4.75e-02 stocproc.method_ft - INFO - J_w_min:1.00e-02 N 64 yields: interval [0.00e+00,4.10e+00] diff 7.88e-03 stocproc.method_ft - INFO - J_w_min:1.00e-04 N 32 yields: interval [0.00e+00,7.50e+00] diff 8.51e-02 stocproc.method_ft - INFO - J_w_min:1.00e-03 N 64 yields: interval [0.00e+00,5.82e+00] diff 1.01e-02 stocproc.method_ft - INFO - J_w_min:1.00e-02 N 128 yields: interval [0.00e+00,4.10e+00] diff 7.56e-03 stocproc.method_ft - INFO - J_w_min:1.00e-05 N 32 yields: interval [0.00e+00,9.16e+00] diff 1.04e-01 stocproc.method_ft - INFO - J_w_min:1.00e-04 N 64 yields: interval [0.00e+00,7.50e+00] diff 1.78e-02 stocproc.method_ft - INFO - J_w_min:1.00e-03 N 128 yields: interval [0.00e+00,5.82e+00] diff 2.38e-03 stocproc.method_ft - INFO - J_w_min:1.00e-02 N 256 yields: interval [0.00e+00,4.10e+00] diff 7.48e-03 stocproc.method_ft - INFO - increasing N while shrinking the interval does lower the error -> try next level stocproc.method_ft - INFO - J_w_min:1.00e-06 N 32 yields: interval [0.00e+00,1.08e+01] diff 1.22e-01 stocproc.method_ft - INFO - J_w_min:1.00e-05 N 64 yields: interval [0.00e+00,9.16e+00] diff 2.81e-02 stocproc.method_ft - INFO - J_w_min:1.00e-04 N 128 yields: interval [0.00e+00,7.50e+00] diff 4.17e-03 stocproc.method_ft - INFO - J_w_min:1.00e-03 N 256 yields: interval [0.00e+00,5.82e+00] diff 7.86e-04 stocproc.method_ft - INFO - return, cause tol of 0.001 was reached stocproc.method_ft - INFO - requires dx < 2.272e-02 stocproc.stocproc - INFO - Fourier Integral Boundaries: [0.000e+00, 7.064e+01] stocproc.stocproc - INFO - Number of Nodes : 4096 stocproc.stocproc - INFO - yields dx : 1.725e-02 stocproc.stocproc - INFO - yields dt : 8.895e-02 stocproc.stocproc - INFO - yields t_max : 3.642e+02 #+end_example The sample trajectories are smooth too. #+begin_src jupyter-python %%space plot ts = np.linspace(0, params.t_max, 1000) EtaTherm.new_process() plt.plot(ts, EtaTherm(ts).real) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/6d082076097694cb45a42d4608ed411592bfcd3e.png]] :END: * Actual Hops Generate the key for binary caching. #+begin_src jupyter-python hi_key = hierarchyData.HIMetaKey_type( HiP=HierachyParam, IntP=IntegrationParam, SysP=SystemParam, Eta=Eta, EtaTherm=EtaTherm, ) #+end_src #+RESULTS: Initialize Hierarchy. #+begin_src jupyter-python myHierarchy = hierarchyLib.HI(hi_key, number_of_samples=params.N, desc="calculate the heat flow") #+end_src #+RESULTS: : init Hi class, use 2002 equation : /home/hiro/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyLib.py:1058: UserWarning: sum_k_max is not implemented! DO SO BEFORE NEXT USAGE (use simplex).HierarchyParametersType does not yet know about sum_k_max : warnings.warn( Run the integration. #+begin_src jupyter-python myHierarchy.integrate_simple(data_path="data", data_name="energy_flow_therm_new_again_less.data") #+end_src #+RESULTS: : samples :[TET 5.48ms [0.0c/s] TTG -- 0.0% ETA -- ORT --] : integration :[TET 5.28ms [0.0c/s] TTG -- 0.0% ETA -- ORT --] : samples :[TET-2.01s---[866.5c/s]-TTG-3.00s-------------> 43.5% ETA 20211103_10:49:30 ORT 5.01s] : integration :[TET 1.71ms [0.0c/s] TTG -- 0.0% ETA -- ORT --] : samples :[TET-4.00s---[999.5c/s]-TTG-0.00ms------------------100%-------------------ETA-20211103_10:49:29-ORT-4.00s] : integration :[TET 3.71ms [0.0c/s] TTG -- 0.0% ETA -- ORT --] :  Get the samples. #+BEGIN_SRC jupyter-python # to access the data the 'hi_key' is used to find the data in the hdf5 file class int_result: with hierarchyData.HIMetaData( hid_name="energy_flow_therm_new_again_less.data", hid_path="data" ) as metaData: with metaData.get_HIData(hi_key, read_only=True) as data: smp = data.get_samples() print("{} samples found in database".format(smp)) τ = data.get_time() rho_τ = data.get_rho_t() #s_proc = np.array(data.stoc_proc) #states = np.array(data.aux_states).copy() ψ_1 = np.array(data.aux_states[:, :, 0 : params.num_exp_t * params.dim]) ψ_0 = np.array(data.stoc_traj) y = np.array(data.y) #η = np.array(data.stoc_proc) temp_y = np.array(data.temp_y) #+end_src #+RESULTS: : 4000 samples found in database Calculate energy. #+begin_src jupyter-python %matplotlib inline energy = np.einsum("ijk,kj", int_result.rho_τ,params.H_s).real plt.plot(int_result.τ, energy) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/8e6031964f047808246e7cc787bff3a3133b953e.png]] :END: * Energy Flow :PROPERTIES: :ID: 9ce93da8-d323-40ec-96a2-42ba184dc963 :END: #+begin_src jupyter-python int_result.ψ_1.shape #+end_src #+RESULTS: | 5120 | 300 | 8 | Let's look at the norm. #+begin_src jupyter-python plt.plot(int_result.τ, (int_result.ψ_0[0].conj() * int_result.ψ_0[0]).sum(axis=1).real) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/46197d5b229e031e905aea1a0f5b517e977f711d.png]] :END: And try to calculate the energy flow. #+begin_src jupyter-python def flow_for_traj(ψ_0, ψ_1, temp_y): a = np.array((params.L @ ψ_0.T).T) EtaTherm.new_process(temp_y) η_dot = scipy.misc.derivative(EtaTherm, int_result.τ, dx=1e-3, order=5) ψ_1 = (-w * g * params.bcf_scale)[None, :, None] * ψ_1.reshape( params.t_steps, params.num_exp_t, params.dim ) # return np.array(np.sum(ψ_0.conj() * ψ_0, axis=1)).flatten().real j_0 = np.array( 2 ,* ( 1j ,* (np.sum(a.conj()[:, None, :] * ψ_1, axis=(1, 2))) / np.sum(ψ_0.conj() * ψ_0, axis=1) ).real ).flatten() j_therm = -np.array( 2 ,* ( (np.sum(a.conj() * ψ_0, axis=1)) * η_dot / np.sum(ψ_0.conj() * ψ_0, axis=1) ).real ).flatten() return j_0, j_therm #+end_src #+RESULTS: Now we calculate the average over all trajectories. #+begin_src jupyter-python class Flow: j_0 = np.zeros_like(int_result.τ) j_therm = np.zeros_like(int_result.τ) for i in range(0, params.N): dj, dj_therm = flow_for_traj( int_result.ψ_0[i], int_result.ψ_1[i], int_result.temp_y[i] ) j_0 += dj j_therm += dj_therm j_0 /= params.N j_therm /= params.N j = j_0 + j_therm #+end_src #+RESULTS: And plot it :). #+begin_src jupyter-python %matplotlib inline plt.plot(int_result.τ, Flow.j_0, label=r"$j_0$") plt.plot(int_result.τ, Flow.j_therm, label=r"$j_\mathrm{therm}$") plt.plot(int_result.τ, Flow.j, label=r"$j$") plt.legend() #+end_src #+RESULTS: :RESULTS: : [[file:./.ob-jupyter/8ee454c5d5d79c6ebc664d4e6c9fa86ddd22d76b.png]] :END: Let's calculate the integrated energy. #+begin_src jupyter-python E_t = np.array( [0] + [ scipy.integrate.simpson(Flow.j[0:n], int_result.τ[0:n]) for n in range(1, len(int_result.τ)) ] ) E_t[-1] #+end_src #+RESULTS: : 0.2055296100424721 With this we can retrieve the energy of the interaction Hamiltonian. #+begin_src jupyter-python E_I = - energy - E_t #+end_src #+RESULTS: #+begin_src jupyter-python %%space plot #plt.plot(τ, j, label="$J$", linestyle='--') plt.plot(int_result.τ, E_t, label=r"$\langle H_{\mathrm{B}}\rangle$") plt.plot(int_result.τ, E_I, label=r"$\langle H_{\mathrm{I}}\rangle$") plt.plot(int_result.τ, energy, label=r"$\langle H_{\mathrm{S}}\rangle$") plt.xlabel("τ") plt.legend() plt.show() #+end_src #+RESULTS: :RESULTS: | | | | | | : Text(0.5, 0, 'τ') : [[file:./.ob-jupyter/4667ed4d0c34454379b42baac0449563ea26df93.png]] :END: * System + Interaction Energy #+begin_src jupyter-python def h_si_for_traj(ψ_0, ψ_1, temp_y): a = np.array((params.L @ ψ_0.T).T) b = np.array((params.H_s @ ψ_0.T).T) ψ_1 = (g*params.bcf_scale)[None, :, None] * ψ_1.reshape( params.t_steps, params.num_exp_t, params.dim ) EtaTherm.new_process(temp_y) E_i = np.array( 2 ,* ( -1j ,* np.sum( a.conj()[:, None, :] ,* ψ_1, axis=(1, 2), ) ).real ).flatten() E_i += np.array( 2 ,* ( EtaTherm(int_result.τ) ,* np.sum( a.conj() ,* ψ_0, axis=1, ) ).real ).flatten() E_s = np.array(np.sum(b.conj() * ψ_0, axis=1)).flatten().real return (E_i + E_s) / np.sum(ψ_0.conj() * ψ_0, axis=1).real #+end_src #+RESULTS: #+begin_src jupyter-python e_si = np.zeros_like(int_result.τ) for i in range(0, params.N): e_si += h_si_for_traj(int_result.ψ_0[i], int_result.ψ_1[i], int_result.temp_y[i]) e_si /= params.N #+end_src #+RESULTS: Doesn't work out. #+begin_src jupyter-python plt.plot(int_result.τ, e_si -energy, label=r"direct") plt.plot(int_result.τ, E_I) plt.legend() #+end_src #+RESULTS: :RESULTS: : [[file:./.ob-jupyter/c8c68dc3ee919f1306e4df7b284c3de1d8d24e8e.png]] :END: