master-thesis/python/richard_hops/energy_flow_nonlinear.org

#+PROPERTY: header-args :session rich_hops_eflow_nl :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
  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
  γ = 3 # coupling ratio
  ω_c = 2  # center of spect. dens
  δ = 2  # breadth BCF
  t_max = 10
  t_steps = 4000
  k_max = 10
  seed = 100
  H_s = σ3 + np.eye(2)
  L =  σ2 * γ # 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

  α = CauchyBCF(δ, ω_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:
| <matplotlib.lines.Line2D | at | 0x7fe0a859b520> |
| <matplotlib.lines.Line2D | at | 0x7fe0a85a88e0> |
| <matplotlib.lines.Line2D | at | 0x7fe0a85a8df0> |
[[file:./.ob-jupyter/e860e1804bf739ce1e809a9560da4176f1f62b3b.png]]
:END:

** Hops setup
#+begin_src jupyter-python
  HierachyParam = hierarchyData.HiP(
      k_max=k_max,
      # g_scale=None,
      # sample_method='random',
      seed=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=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_FFT(
      α.I,
      t_max,
      α,
      negative_frequencies=True,
      seed=seed,
      intgr_tol=1e-2,
      intpl_tol=1e-2,
      scale=1,
  )
#+end_src

#+RESULTS:
#+begin_example
  stocproc.stocproc - INFO - use neg freq
  stocproc.method_ft - INFO - get_dt_for_accurate_interpolation, please wait ...
  stocproc.method_ft - INFO - acc interp N 33 dt 1.44e-01 -> diff 7.57e-03
  stocproc.method_ft - INFO - requires dt < 1.439e-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 [-1.47e+01,1.87e+01] diff 3.37e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 32 yields: interval [-5.11e+01,5.51e+01] diff 6.70e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-02 N 64 yields: interval [-1.47e+01,1.87e+01] diff 3.37e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 32 yields: interval [-1.66e+02,1.70e+02] diff 2.44e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 64 yields: interval [-5.11e+01,5.51e+01] diff 1.11e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-02 N 128 yields: interval [-1.47e+01,1.87e+01] diff 3.37e-01
  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 [-5.30e+02,5.34e+02] diff 3.68e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 64 yields: interval [-1.66e+02,1.70e+02] diff 1.34e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 128 yields: interval [-5.11e+01,5.51e+01] diff 1.06e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-02 N 256 yields: interval [-1.47e+01,1.87e+01] diff 3.37e-01
  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 [-1.68e+03,1.68e+03] diff 4.19e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-05 N 64 yields: interval [-5.30e+02,5.34e+02] diff 3.04e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 128 yields: interval [-1.66e+02,1.70e+02] diff 4.07e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 256 yields: interval [-5.11e+01,5.51e+01] diff 1.06e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-02 N 512 yields: interval [-1.47e+01,1.87e+01] diff 3.37e-01
  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-07 N 32 yields: interval [-5.32e+03,5.32e+03] diff 4.36e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-06 N 64 yields: interval [-1.68e+03,1.68e+03] diff 3.94e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-05 N 128 yields: interval [-5.30e+02,5.34e+02] diff 2.09e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 256 yields: interval [-1.66e+02,1.70e+02] diff 3.72e-02
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 512 yields: interval [-5.11e+01,5.51e+01] diff 1.06e-01
  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-08 N 32 yields: interval [-1.68e+04,1.68e+04] diff 4.42e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-07 N 64 yields: interval [-5.32e+03,5.32e+03] diff 4.28e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-06 N 128 yields: interval [-1.68e+03,1.68e+03] diff 3.50e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-05 N 256 yields: interval [-5.30e+02,5.34e+02] diff 9.79e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 512 yields: interval [-1.66e+02,1.70e+02] diff 3.36e-02
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 1024 yields: interval [-5.11e+01,5.51e+01] diff 1.06e-01
  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-09 N 32 yields: interval [-5.32e+04,5.32e+04] diff 4.43e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-08 N 64 yields: interval [-1.68e+04,1.68e+04] diff 4.39e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-07 N 128 yields: interval [-5.32e+03,5.32e+03] diff 4.12e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-06 N 256 yields: interval [-1.68e+03,1.68e+03] diff 2.75e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-05 N 512 yields: interval [-5.30e+02,5.34e+02] diff 2.16e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 1024 yields: interval [-1.66e+02,1.70e+02] diff 3.36e-02
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 2048 yields: interval [-5.11e+01,5.51e+01] diff 1.06e-01
  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-10 N 32 yields: interval [-1.68e+05,1.68e+05] diff 4.44e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-09 N 64 yields: interval [-5.32e+04,5.32e+04] diff 4.43e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-08 N 128 yields: interval [-1.68e+04,1.68e+04] diff 4.34e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-07 N 256 yields: interval [-5.32e+03,5.32e+03] diff 3.82e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-06 N 512 yields: interval [-1.68e+03,1.68e+03] diff 1.71e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-05 N 1024 yields: interval [-5.30e+02,5.34e+02] diff 1.07e-02
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 2048 yields: interval [-1.66e+02,1.70e+02] diff 3.36e-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-11 N 32 yields: interval [-5.32e+05,5.32e+05] diff 4.44e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-10 N 64 yields: interval [-1.68e+05,1.68e+05] diff 4.44e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-09 N 128 yields: interval [-5.32e+04,5.32e+04] diff 4.41e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-08 N 256 yields: interval [-1.68e+04,1.68e+04] diff 4.24e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-07 N 512 yields: interval [-5.32e+03,5.32e+03] diff 3.28e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-06 N 1024 yields: interval [-1.68e+03,1.68e+03] diff 6.56e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-05 N 2048 yields: interval [-5.30e+02,5.34e+02] diff 1.06e-02
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 4096 yields: interval [-1.66e+02,1.70e+02] diff 3.36e-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-12 N 32 yields: interval [-1.68e+06,1.68e+06] diff 4.44e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-11 N 64 yields: interval [-5.32e+05,5.32e+05] diff 4.44e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-10 N 128 yields: interval [-1.68e+05,1.68e+05] diff 4.43e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-09 N 256 yields: interval [-5.32e+04,5.32e+04] diff 4.38e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-08 N 512 yields: interval [-1.68e+04,1.68e+04] diff 4.04e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-07 N 1024 yields: interval [-5.32e+03,5.32e+03] diff 2.43e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-06 N 2048 yields: interval [-1.68e+03,1.68e+03] diff 9.69e-02
  stocproc.method_ft - INFO - J_w_min:1.00e-05 N 4096 yields: interval [-5.30e+02,5.34e+02] diff 1.06e-02
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 8192 yields: interval [-1.66e+02,1.70e+02] diff 3.36e-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-13 N 32 yields: interval [-5.32e+06,5.32e+06] diff 4.44e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-12 N 64 yields: interval [-1.68e+06,1.68e+06] diff 4.44e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-11 N 128 yields: interval [-5.32e+05,5.32e+05] diff 4.44e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-10 N 256 yields: interval [-1.68e+05,1.68e+05] diff 4.42e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-09 N 512 yields: interval [-5.32e+04,5.32e+04] diff 4.31e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-08 N 1024 yields: interval [-1.68e+04,1.68e+04] diff 3.67e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-07 N 2048 yields: interval [-5.32e+03,5.32e+03] diff 1.33e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-06 N 4096 yields: interval [-1.68e+03,1.68e+03] diff 3.37e-03
  stocproc.method_ft - INFO - return, cause tol of 0.01 was reached
  stocproc.method_ft - INFO - requires dx < 8.212e-01
  stocproc.method_ft - INFO - increase N to match dt_new*(N-1) < t_max
  stocproc.stocproc - INFO - Fourier Integral Boundaries: [-1.680e+03, 1.684e+03]
  stocproc.stocproc - INFO - Number of Nodes            : 8192
  stocproc.stocproc - INFO - yields dx                  : 4.106e-01
  stocproc.stocproc - INFO - yields dt                  : 1.868e-03
  stocproc.stocproc - INFO - yields t_max               : 1.530e+01
#+end_example

* 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 22 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_name="energy_flow_nl.data")
#+end_src

#+RESULTS:
:RESULTS:
#+begin_example
  (10, 4000, 11)
  samples     :0.0%
  integration :0.0%
  samples     :0.0%
  integration :22.0%
  samples     :0.0%
  integration :43.7%
  samples     :0.0%
  integration :65.7%
  samples     :0.0%
  integration :87.6%
  samples     :1.0%
  integration :11.5%
  Process Process-591:
  Traceback (most recent call last):
  File "/nix/store/75svh7wjqcdc21x23xm99zq8fckqmgxb-python3-3.9.4/lib/python3.9/multiprocessing/process.py", line 318, in _bootstrap
    util._exit_function()
  File "/nix/store/75svh7wjqcdc21x23xm99zq8fckqmgxb-python3-3.9.4/lib/python3.9/multiprocessing/util.py", line 334, in _exit_function
    _run_finalizers(0)
  File "/nix/store/75svh7wjqcdc21x23xm99zq8fckqmgxb-python3-3.9.4/lib/python3.9/multiprocessing/util.py", line 291, in _run_finalizers
    keys = [key for key in list(_finalizer_registry) if f(key)]
  File "/home/hiro/Documents/Projects/UNI/master/masterarb/python/richard_hops/progression/progress.py", line 1312, in _stop_on_signal
    raise LoopInterruptError()
  progression.progress.LoopInterruptError
  samples     :1.0%
  integration :25.9%
  
#+end_example
# [goto error]
#+begin_example

  KeyboardInterruptTraceback (most recent call last)
  <ipython-input-913-cc745386a46d> in <module>
  ----> 1 myHierarchy.integrate_simple(data_name="energy_flow_nl.data")

  ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyLib.py in integrate_simple(self, data_name, data_path, overwrite, clear_pd)
     1308
     1309                     if do_calculation:
  -> 1310                         (t, psi_all, e_and_trb), z = _integrate_hierarchy(
     1311                             arg, const_arg, c_int, m_int
     1312                         )

  ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyLib.py in _integrate_hierarchy(arg, const_arg, c, m)
      873
      874     # t0, t1, N, f, args, x0, integrator, verbose, c, **kwargs
  --> 875     return ode_wrapper.integrate_cplx(c=c, args=args_dgl, **kwargs), z
      876
      877

  ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/ode_wrapper.py in integrate_cplx(c, t0, t1, N, f, args, x0, integrator, res_dim, x_to_res, **kwargs)
      123                 while r.successful() and i < N:
      124                     _t = time()
  --> 125                     r.integrate(t[i])
      126                     t_int += (time()-_t)
      127

  /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/scipy/integrate/_ode.py in integrate(self, t, step, relax)
      431
      432         try:
  --> 433             self._y, self.t = mth(self.f, self.jac or (lambda: None),
      434                                   self._y, self.t, t,
      435                                   self.f_params, self.jac_params)

  /nix/store/7r8xg0344zc6lhyqqk2lynwbh8hy3934-python3-3.9.4-env/lib/python3.9/site-packages/scipy/integrate/_ode.py in run(self, f, jac, y0, t0, t1, f_params, jac_params)
     1007         args = ((f, jac, y0, t0, t1) + tuple(self.call_args) +
     1008                 (f_params, jac_params))
  -> 1009         y1, t, istate = self.runner(*args)
     1010         self.istate = istate
     1011         if istate < 0:

  ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/ode_wrapper.py in new_f(t, x)
       48     def new_f(t, x):
       49         t0 = time()
  ---> 50         res = f(t, x)
       51         t1 = time()
       52         time_as_list[0] += t1-t0

  ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/ode_wrapper.py in <lambda>(t, x)
       56
       57 def integrate_cplx(c, t0, t1, N, f, args, x0, integrator, res_dim=None, x_to_res=None, **kwargs):
  ---> 58     f_partial_complex = lambda t, x: f(t, x, *args)
       59     if integrator == 'zvode':
       60         # define complex derivative

  ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyLib.py in f_non_lin_hierarchy(t, V_psi_and_eta_lambda, eta_stoc, O_vec, K, B, C, M_down, M_up, spn)
      385
      386     ddt_eta_lambda = ddt_eta_lambda_non_lin(spn, V_psi, eta_lambda)
  --> 387     ddt_psi = f_hierarchy(
      388         t, V_psi, eta_stoc, eta_det, O_vec, K, B, C, M_down, M_up, spn
      389     )

  ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyLib.py in f_hierarchy(t, V_psi, eta_stoc, eta_det, O_vec, K, B, C, M_down, M_up, spn)
      239         O_vec * psi
      240         + K_mat.dot(psi.T).T
  --> 241         + B(t, psi, spn).dot(M_up.dot(psi).T).T
      242         + C(t, psi, spn).dot(M_down.dot(psi).T).T
      243     )

  ~/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyLib.py in B_non_lin(t, psi, spn)
      501     """
      502     psi_zero = psi[0, :]
  --> 503     return spn.minus_L_dagger + E(spn, psi_zero) * spn.eye
      504
      505

  KeyboardInterrupt:
#+end_example
:END:
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_nl.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)
          η = np.array(data.stoc_proc)
#+end_src

#+RESULTS:
: 100 samples found in database

Calculate energy.
#+begin_src jupyter-python
  %matplotlib inline
  energy = np.array([np.trace(ρ @ H_s).real for ρ in rho_τ])
  plt.plot(τ, energy)
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7fe0a8cc82b0> |
[[file:./.ob-jupyter/fc86a420fa0dde236981598fd9e3332741a55ab5.png]]
:END:

#+begin_src jupyter-python
  %%space plot
  plt.plot(τ, np.trace(rho_τ.T).real)
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7fe0a8c9dc40> |
[[file:./.ob-jupyter/23900676701a8e2760d8e6d0d758bc8b2d27f7f7.png]]
:END:

* Energy Flow
:PROPERTIES:
:ID:       eefb1594-e399-4d24-9dd7-a57addd42e65
:END:
#+begin_src jupyter-python
  ψ_1.shape
#+end_src

#+RESULTS:
| 160 | 4000 | 2 |

Let's look at the norm.
#+begin_src jupyter-python
  plt.plot(τ, (ψ_0[0].conj() * ψ_0[0]).sum(axis=1).real)
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7fe0a8c0e790> |
[[file:./.ob-jupyter/122288a6cbed863e55fe47b9cd4590bf5bcd6207.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(np.sum(ψ_0.conj() * ψ_0, axis=1)).flatten().real
      return np.array(2 * (1j * -W * np.sum(a.conj() * ψ_1, axis=1)/np.sum(ψ_0.conj() * ψ_0, 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:

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:

And plot it :)
#+begin_src jupyter-python
  %matplotlib inline
  plt.plot(τ, j)
  #plt.plot(τ, ja)
  plt.show()
#+end_src

#+RESULTS:
[[file:./.ob-jupyter/2ced0198134dfd75dbcd914579d64b75fba0ab86.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:
: 4.559650096824831

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'] = [15, 10]
  #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:
| <matplotlib.lines.Line2D | at | 0x7fe0a8b7e1c0> |
| <matplotlib.lines.Line2D | at | 0x7fe0a8b7e610> |
| <matplotlib.lines.Line2D | at | 0x7fe0a8b7e9a0> |
: Text(0.5, 0, 'τ')
: <matplotlib.legend.Legend at 0x7fe0a8b7ebb0>
[[file:./.ob-jupyter/f1f59306d9a3a957e251ba42e32ee5907656c1a5.png]]
:END:
#+RESULTS:

* System + Interaction Energy
#+begin_src jupyter-python
  def h_si_for_traj(ψ_0, ψ_1):
      a = np.array((L @ ψ_0.T).T)
      b = np.array((H_s @ ψ_0.T).T)
      E_i = np.array(2 * (-1j * np.sum(a.conj() * ψ_1, 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

  def h_si_for_traj_alt(ψ_0, y):
      Eta.new_process(y)

      a = np.array((L.conj().T @ ψ_0.T).T)
      b = np.array((H_s @ ψ_0.T).T)
      E_i = np.array(2 * (Eta(τ) * 1j * 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
#+end_src

#+RESULTS:

#+begin_src jupyter-python
  e_si = np.zeros_like(τ)
  for i in range(0, N):
      e_si += h_si_for_traj(ψ_0[i], ψ_1[i])
  e_si /= N
#+end_src

#+RESULTS:

Not too bad...
#+begin_src jupyter-python
  plt.plot(τ, e_si)
  plt.ylim(-16,2)
  plt.plot(τ, E_I + energy)
#+end_src

#+RESULTS:
:RESULTS:
| <matplotlib.lines.Line2D | at | 0x7fe0a866afa0> |
[[file:./.ob-jupyter/0ce4e1fbfe77b4a27d7b2f7ec1978a99e6974805.png]]
:END: