#+PROPERTY: header-args :session rich_hops_eflow :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 : The jupyter_spaces extension is already loaded. To reload it, use: : %reload_ext jupyter_spaces ** 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, StocProc_KLE import bcf from dataclasses import dataclass import scipy import scipy.misc import scipy.signal #+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 γ = 5 # coupling ratio ω_c = 0 # center of spect. dens δ = .1 # breadth BCF t_max = 10 t_steps = 500 k_max = 6 seed = 100 H_s = σ3 + np.eye(2) L = 1 / 2 * (σ1 - 1j * σ2) * γ ψ_0 = np.array([1, 0]) W = ω_c * 1j + δ # exponent BCF N = 100 #+end_src #+RESULTS: ** BCF #+begin_src jupyter-python @dataclass class CauchyBCF: δ: float ω_c: float def I(self, ω): return np.sqrt(self.δ) / (self.δ + (ω - self.ω_c) ** 2 / self.δ) def __call__(self, τ): return np.sqrt(self.δ) * np.exp(-1j * self.ω_c * τ - np.abs(τ) * self.δ) def __bfkey__(self): return self.δ, self.ω_c @dataclass class GaussBCF: σ: float ω_c: float def I(self, ω): return ( np.exp(-(((ω - self.ω_c) / self.σ) ** 2) / 2) ,* 1 / (np.sqrt(2 * np.pi) * self.σ) * np.pi ) def __call__(self, τ): return ( np.exp(-(((τ - self.ω_c) / self.σ) ** 2) / 2) ,* 1 / (np.sqrt(2 * np.pi) * self.σ) * np.pi ) #np.exp(1j * self.ω_c * τ - self.σ**2 * τ**2/2) def __bfkey__(self): return self.σ, self.ω_c α = GaussBCF(δ, ω_c) #+end_src #+RESULTS: *** Plot #+begin_src jupyter-python %%space plot t = np.linspace(0, t_max, 1000) ω = np.linspace(ω_c - 10, ω_c + 10, 1000) fig, axs = plt.subplots(2) axs[0].plot(t, np.real(α(t))) axs[0].plot(t, np.imag(α(t))) axs[1].plot(ω, α.I(ω)) #+end_src #+RESULTS: :RESULTS: | | | | | | [[file:./.ob-jupyter/e79e0652b99bba8df94a7f60c64af072947ade03.png]] :END: ** Hops setup #+begin_src jupyter-python HierachyParam = hierarchyData.HiP( k_max=k_max, # g_scale=None, # sample_method='random', seed=seed, nonlinear=False, # 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=t_max, t_steps=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=H_s, L=L, psi0=ψ_0, # excited qubit g=np.array([np.sqrt(δ)]), w=np.array([W]), H_dynamic=[], bcf_scale=1, # 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=1, ) #+end_src #+RESULTS: The quantum noise. #+begin_src jupyter-python Eta = StocProc_KLE( α, t_max, seed=seed, tol=1e-3 ) #+end_src #+RESULTS: :RESULTS: #+begin_example stocproc.method_kle - INFO - check 33 grid points stocproc.method_kle - INFO - check 65 grid points alpha_k is real alpha_k is real stocproc.method_kle - INFO - check 129 grid points alpha_k is real stocproc.method_kle - INFO - check 257 grid points alpha_k is real stocproc.method_kle - INFO - check 513 grid points alpha_k is real #+end_example # [goto error] #+begin_example KeyboardInterruptTraceback (most recent call last) in ----> 1 Eta = StocProc_KLE( 2 α, 3 t_max, 4 seed=seed, 5 tol=1e-3 /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/stocproc/stocproc.py in __init__(self, alpha, t_max, tol, ng_fac, meth, diff_method, dm_random_samples, seed, align_eig_vec, scale) 286 key = alpha, t_max, tol 287 --> 288 sqrt_lambda_ui_fine, t = method_kle.auto_ng( 289 corr=alpha, 290 t_max=t_max, /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/stocproc/method_kle.py in auto_ng(corr, t_max, ngfac, meth, tol, diff_method, dm_random_samples, ret_eigvals, relative_difference) 543 ui_super_fine.reshape(1, -1) 544 ) --> 545 md = np.max(np.abs(diff) / abs_alpha_res) 546 time_calc_diff += time.time() - t0 547 KeyboardInterrupt: #+end_example :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=None, ) #+end_src #+RESULTS: Initialize Hierarchy. #+begin_src jupyter-python myHierarchy = hierarchyLib.HI(hi_key, number_of_samples=N, desc="run a test case") #+end_src #+RESULTS: : init Hi class, use 14 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 :results none myHierarchy.integrate_simple(data_name="energy_flow.data") #+end_src Get the samples. #+begin_src jupyter-python # to access the data the 'hi_key' is used to find the data in the hdf5 file with hierarchyData.HIMetaData(hid_name="energy_flow.data", hid_path=".") 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:2] ψ_0 = np.array(data.stoc_traj) y = np.array(data.y) #+end_src #+RESULTS: :RESULTS: # [goto error] #+begin_example KeyErrorTraceback (most recent call last) ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyData.py in get_HIData_fname(self, key) 779 try: --> 780 hdf5_name = self.db[hashed_key][0] 781 except KeyError: /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/sqlitedict.py in __getitem__(self, key) 243 if item is None: --> 244 raise KeyError(key) 245 return self.decode(item[0]) KeyError: '47177f42579f772ba59f8489393910e4b2e9b1f2567e082f6b944d00382a9df7793b33c105854cabecd012bc9f3fb59d' During handling of the above exception, another exception occurred: PicklingErrorTraceback (most recent call last) in 1 # to access the data the 'hi_key' is used to find the data in the hdf5 file 2 with hierarchyData.HIMetaData(hid_name="energy_flow.data", hid_path=".") as metaData: ----> 3 with metaData.get_HIData(hi_key, read_only=True) as data: 4 smp = data.get_samples() 5 print("{} samples found in database".format(smp)) ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyData.py in get_HIData(self, key, read_only) 787 788 def get_HIData(self, key, read_only=False): --> 789 hdf5_name = self.get_HIData_fname(key) 790 791 if key.HiP.result_type == RESULT_TYPE_ZEROTH_ORDER_ONLY: ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyData.py in get_HIData_fname(self, key) 781 except KeyError: 782 hdf5_name = self._new_rand_file_name(pre=self.name + "_", end=".h5") --> 783 self.db[hashed_key] = (hdf5_name, key) 784 self.db.commit() 785 /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/sqlitedict.py in __setitem__(self, key, value) 250 251 ADD_ITEM = 'REPLACE INTO "%s" (key, value) VALUES (?,?)' % self.tablename --> 252 self.conn.execute(ADD_ITEM, (key, self.encode(value))) 253 if self.autocommit: 254 self.commit() /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/sqlitedict.py in encode(obj) 95 def encode(obj): 96 """Serialize an object using pickle to a binary format accepted by SQLite.""" ---> 97 return sqlite3.Binary(dumps(obj, protocol=PICKLE_PROTOCOL)) 98 99 PicklingError: Can't pickle : it's not the same object as __main__.GaussBCF #+end_example :END: Calculate energy. #+begin_src jupyter-python energy = np.array([np.trace(ρ * H_s).real/np.trace(ρ).real for ρ in rho_τ]) plt.plot(τ, energy) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/892241bc3127d1dc4581b7f31e8b392425c9716c.png]] :END: #+begin_src jupyter-python %%space plot plt.plot(τ, np.trace(rho_τ.T).real) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/d79cca7773c99c7279d79f331d77f301acbdb71d.png]] :END: * Energy Flow :PROPERTIES: :ID: eefb1594-e399-4d24-9dd7-a57addd42e65 :END: #+begin_src jupyter-python ψ_1.shape #+end_src #+RESULTS: | 160 | 500 | 2 | Let's look at the norm. #+begin_src jupyter-python plt.plot(τ, (ψ_1[0].conj() * ψ_1[0]).sum(axis=1).real) #+end_src #+RESULTS: :RESULTS: | | [[file:./.ob-jupyter/5783eff51f0815225fe39255d41152e4b447963f.png]] :END: And try to calculate the energy flow. #+begin_src jupyter-python def flow_for_traj(ψ_0, ψ_1): a = np.array((L @ ψ_0.T).T) return np.array(2 * (1j * -W * np.sum(a.conj() * ψ_1, axis=1)).real).flatten() def flow_for_traj_alt(ψ_0, y): Eta.new_process(y) eta_dot = scipy.misc.derivative(Eta, τ, dx=1e-8) a = np.array((L @ ψ_0.T).T) return -(2j * eta_dot.conj() * np.array((np.sum(ψ_0.conj() * a, axis=1))).flatten() ).real #+end_src #+RESULTS: Now we calculate the average over all trajectories. #+begin_src jupyter-python j = np.zeros_like(τ) for i in range(0, N): j += flow_for_traj(ψ_0[i], ψ_1[i]) j /= N #+end_src #+RESULTS: :RESULTS: # [goto error] #+begin_example ValueErrorTraceback (most recent call last) in 1 j = np.zeros_like(τ) 2 for i in range(0, N): ----> 3 j += flow_for_traj(ψ_0[i], ψ_1[i]) 4 j /= N in flow_for_traj(ψ_0, ψ_1) 1 def flow_for_traj(ψ_0, ψ_1): ----> 2 a = np.array((L @ ψ_0.T).T) 3 4 return np.array(2 * (1j * -W * np.sum(a.conj() * ψ_1, axis=1)).real).flatten() 5 ValueError: matmul: Input operand 1 does not have enough dimensions (has 0, gufunc core with signature (n?,k),(k,m?)->(n?,m?) requires 1) #+end_example :END: And do the same with the alternative implementation. #+begin_src jupyter-python ja = np.zeros_like(τ) for i in range(0, N): ja += flow_for_traj_alt(ψ_0[i], y[i]) ja /= N #+end_src #+RESULTS: :RESULTS: # [goto error] #+begin_example ValueErrorTraceback (most recent call last) in 1 ja = np.zeros_like(τ) 2 for i in range(0, N): ----> 3 ja += flow_for_traj_alt(ψ_0[i], y[i]) 4 ja /= N in flow_for_traj_alt(ψ_0, y) 8 Eta.new_process(y) 9 eta_dot = scipy.misc.derivative(Eta, τ, dx=1e-8) ---> 10 a = np.array((L @ ψ_0.T).T) 11 12 return -(2j * eta_dot.conj() * ValueError: matmul: Input operand 1 does not have enough dimensions (has 0, gufunc core with signature (n?,k),(k,m?)->(n?,m?) requires 1) #+end_example :END: And plot it :) #+begin_src jupyter-python %matplotlib inline plt.plot(τ, j) plt.plot(τ, ja) plt.show() #+end_src #+RESULTS: [[file:./.ob-jupyter/b159f74a1bffc1460ac7f285f76ab1b81f31bd07.png]] \Let's calculate the integrated energy. #+begin_src jupyter-python E_t = np.array([0] + [scipy.integrate.simpson(j[0:n], τ[0:n]) for n in range(1, len(τ))]) E_t[-1] #+end_src #+RESULTS: : 0.0 #+begin_src jupyter-python E_t = np.array([0] + [scipy.integrate.simpson(ja[0:n], τ[0:n]) for n in range(1, len(τ))]) E_t[-1] #+end_src #+RESULTS: : 0.0 With this we can retrieve the energy of the interaction Hamiltonian. #+begin_src jupyter-python E_I = 2 - energy - E_t #+end_src #+RESULTS: #+begin_src jupyter-python %%space plot plt.rcParams['figure.figsize'] = [10, 8] #plt.plot(τ, j, label="$J$", linestyle='--') plt.plot(τ, E_t, label=r"$\langle H_{\mathrm{B}}\rangle$") plt.plot(τ, E_I, label=r"$\langle H_{\mathrm{I}}\rangle$") plt.plot(τ, 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/ca75f76563aca062310d5779c0c5df539b728d3f.png]] :END: #+RESULTS: * Derivatives #+begin_src jupyter-python Eta.new_process(y[0]) #plt.plot(τ, Eta(τ).real) eta_dot = scipy.misc.derivative(Eta, τ, dx=1e-3) plt.plot(τ, eta_dot) #+end_src #+RESULTS: :RESULTS: : /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/numpy/core/_asarray.py:102: ComplexWarning: Casting complex values to real discards the imaginary part : return array(a, dtype, copy=False, order=order) | | [[file:./.ob-jupyter/aaa6bdcc61452758e092eccc2d62ef8cdfe5f258.png]] :END: