master-thesis/python/richard_hops/energy_flow.org

• Setup
  • Jupyter
  • Matplotlib
  • Richard (old) HOPS
  • Auxiliary Definitions
• Model Setup
  • BCF
    • Plot
  • Hops setup
• Actual Hops
• Energy Flow

Setup

Jupyter

  %load_ext autoreload
  %autoreload 2
  %load_ext jupyter_spaces

Matplotlib

  import matplotlib
  import matplotlib.pyplot as plt

  #matplotlib.use("TkCairo", force=True)
  %gui tk
  %matplotlib inline
  plt.style.use('ggplot')

Richard (old) HOPS

  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

Auxiliary Definitions

  σ1 = np.matrix([[0,1],[1,0]])
  σ2 = np.matrix([[0,-1j],[1j,0]])
  σ3 = np.matrix([[1,0],[0,-1]])

Model Setup

Basic parameters.

  γ = 5 # coupling ratio
  ω_c = 1  # center of spect. dens
  δ = 1  # breadth BCF
  t_max = 10
  t_steps = 500
  k_max = 2
  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 = 10

BCF

  @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)

Plot

  %%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(ω))

<matplotlib.lines.Line2D	at	0x7efc7a2362e0>
<matplotlib.lines.Line2D	at	0x7efc7a236760>
<matplotlib.lines.Line2D	at	0x7efc7a236b50>

Hops setup

  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
  )

Integration.

  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
  )

And now the system.

  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,
  )

The quantum noise.

  Eta = StocProc_FFT(
      α.I,
      t_max,
      α,
      negative_frequencies=True,
      seed=seed,
      intgr_tol=1e-2,
      intpl_tol=1e-2,
      scale=1,
  )

  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 2.88e-01 -> diff 7.57e-03
  stocproc.method_ft - INFO - requires dt < 2.878e-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 [-8.95e+00,1.09e+01] diff 2.01e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 32 yields: interval [-3.06e+01,3.26e+01] diff 6.40e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-02 N 64 yields: interval [-8.95e+00,1.09e+01] diff 2.00e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 32 yields: interval [-9.90e+01,1.01e+02] diff 1.90e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 64 yields: interval [-3.06e+01,3.26e+01] diff 1.31e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-02 N 128 yields: interval [-8.95e+00,1.09e+01] diff 2.00e-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 [-3.15e+02,3.17e+02] diff 2.68e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 64 yields: interval [-9.90e+01,1.01e+02] diff 1.15e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 128 yields: interval [-3.06e+01,3.26e+01] diff 6.33e-02
  stocproc.method_ft - INFO - J_w_min:1.00e-02 N 256 yields: interval [-8.95e+00,1.09e+01] diff 2.00e-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 [-9.99e+02,1.00e+03] diff 2.99e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-05 N 64 yields: interval [-3.15e+02,3.17e+02] diff 2.29e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 128 yields: interval [-9.90e+01,1.01e+02] diff 4.21e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 256 yields: interval [-3.06e+01,3.26e+01] diff 6.33e-02
  stocproc.method_ft - INFO - J_w_min:1.00e-02 N 512 yields: interval [-8.95e+00,1.09e+01] diff 2.00e-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 [-3.16e+03,3.16e+03] diff 3.09e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-06 N 64 yields: interval [-9.99e+02,1.00e+03] diff 2.84e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-05 N 128 yields: interval [-3.15e+02,3.17e+02] diff 1.66e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 256 yields: interval [-9.90e+01,1.01e+02] diff 5.63e-02
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 512 yields: interval [-3.06e+01,3.26e+01] diff 6.33e-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-08 N 32 yields: interval [-1.00e+04,1.00e+04] diff 3.13e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-07 N 64 yields: interval [-3.16e+03,3.16e+03] diff 3.04e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-06 N 128 yields: interval [-9.99e+02,1.00e+03] diff 2.57e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-05 N 256 yields: interval [-3.15e+02,3.17e+02] diff 8.81e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 512 yields: interval [-9.90e+01,1.01e+02] diff 2.00e-02
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 1024 yields: interval [-3.06e+01,3.26e+01] diff 6.33e-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-09 N 32 yields: interval [-3.16e+04,3.16e+04] diff 3.14e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-08 N 64 yields: interval [-1.00e+04,1.00e+04] diff 3.11e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-07 N 128 yields: interval [-3.16e+03,3.16e+03] diff 2.95e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-06 N 256 yields: interval [-9.99e+02,1.00e+03] diff 2.10e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-05 N 512 yields: interval [-3.15e+02,3.17e+02] diff 2.47e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 1024 yields: interval [-9.90e+01,1.01e+02] diff 2.00e-02
  stocproc.method_ft - INFO - J_w_min:1.00e-03 N 2048 yields: interval [-3.06e+01,3.26e+01] diff 6.33e-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-10 N 32 yields: interval [-1.00e+05,1.00e+05] diff 3.14e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-09 N 64 yields: interval [-3.16e+04,3.16e+04] diff 3.13e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-08 N 128 yields: interval [-1.00e+04,1.00e+04] diff 3.08e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-07 N 256 yields: interval [-3.16e+03,3.16e+03] diff 2.77e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-06 N 512 yields: interval [-9.99e+02,1.00e+03] diff 1.41e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-05 N 1024 yields: interval [-3.15e+02,3.17e+02] diff 1.94e-02
  stocproc.method_ft - INFO - J_w_min:1.00e-04 N 2048 yields: interval [-9.90e+01,1.01e+02] diff 2.00e-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 [-3.16e+05,3.16e+05] diff 3.14e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-10 N 64 yields: interval [-1.00e+05,1.00e+05] diff 3.14e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-09 N 128 yields: interval [-3.16e+04,3.16e+04] diff 3.12e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-08 N 256 yields: interval [-1.00e+04,1.00e+04] diff 3.02e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-07 N 512 yields: interval [-3.16e+03,3.16e+03] diff 2.44e+00
  stocproc.method_ft - INFO - J_w_min:1.00e-06 N 1024 yields: interval [-9.99e+02,1.00e+03] diff 6.29e-01
  stocproc.method_ft - INFO - J_w_min:1.00e-05 N 2048 yields: interval [-3.15e+02,3.17e+02] diff 6.32e-03
  stocproc.method_ft - INFO - return, cause tol of 0.01 was reached
  stocproc.method_ft - INFO - requires dx < 3.088e-01
  stocproc.stocproc - INFO - Fourier Integral Boundaries: [-3.152e+02, 3.172e+02]
  stocproc.stocproc - INFO - Number of Nodes            : 2048
  stocproc.stocproc - INFO - yields dx                  : 3.088e-01
  stocproc.stocproc - INFO - yields dt                  : 9.935e-03
  stocproc.stocproc - INFO - yields t_max               : 2.034e+01

Actual Hops

Generate the key for binary caching.

  hi_key = hierarchyData.HIMetaKey_type(
      HiP=HierachyParam,
      IntP=IntegrationParam,
      SysP=SystemParam,
      Eta=Eta,
      EtaTherm=None,
  )

Initialize Hierarchy.

  myHierarchy = hierarchyLib.HI(hi_key, number_of_samples=N, desc="run a test case")

init Hi class, use 6 equation
/home/hiro/Documents/Projects/UNI/master/masterarb/python/richard_hops/hierarchyLib.py:998: 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('sum_k_max is not implemented! DO SO BEFORE NEXT USAGE (use simplex).'+

Run the integration.

  myHierarchy.integrate_simple(data_name="energy_flow.data", overwrite=True)

  samples     :0.0%
  integration :0.0%
  [-1.21674474-0.53065229j  1.29843781-0.42212396j -0.06045702+0.81178203j
   ... -1.03648954-0.03013739j -0.39385685-1.0388868j
   -0.18852159+0.40768442j]
  [ 0.36097777-0.20999212j -0.76016709-0.32886212j  0.20155971+0.78786114j
   ... -0.20699756+0.56375004j -0.21725595+0.47253688j
    0.80215377+0.89727163j]
  samples     :20.0%
  integration :42.0%
  [-1.36534845-0.26838665j  0.10660685-0.46145127j -0.61427558-0.13056982j
  ... -1.20681667-0.79059447j -0.4072043 -0.22122183j
   0.04274245+0.82443994j]
  [-1.22861576-0.13916935j  0.0218835 -0.9336182j   0.01948521+1.03111422j
   ...  0.86619701-0.72840755j -0.41004916+0.5181127j
   -0.05269705+0.82821403j]
  [0.08247194-0.95957482j 0.25563361-0.00802737j 0.45539223-0.70437692j ...
   0.39605661+0.40318642j 0.5208233 -0.74711934j 0.19058665-0.38544974j]
  samples     :50.0%
  integration :17.8%
  [ 0.44834688+0.26218914j -0.95569788+0.00697922j -0.9354395 +0.7900597j
  ...  0.15016663-0.55217627j -0.85241527+0.96090231j
  -0.59842916+0.75357213j]
  [-0.49081064+0.25518693j -0.30969903+0.11565907j  0.86979489-0.67645523j
   ... -1.13665739+0.0503456j   0.0348233 +0.37451812j
    0.86601417-0.16270255j]
  [-1.0160508 -0.16911144j  0.5200195 -1.00714817j -1.05635398-0.30223967j
   ... -0.08468235+0.24238194j  0.38477837-0.40781186j
   -0.53370329-0.72333536j]
  samples     :80.0%
  integration :4.4%
  [-0.42590091+0.75942462j  0.24692167-0.93364911j -1.01259496-0.013466j
  ... -0.28261085+0.11964996j  1.19635194-0.16821297j
  -0.06671849+1.11275477j]
  [ 0.1282879 +0.52498879j  0.28174411-0.5695823j  -0.41456977-0.48066038j
   ...  0.51724726-0.72920012j  0.76341678+1.14227298j
    0.02779715-0.1876853j ]
  samples     : 100%
  integration :0.0%
  

Get the samples.

  # 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)

122 samples found in database

Calculate energy.

  energy = np.array([np.trace(ρ * H_s).real/np.trace(ρ).real for ρ in rho_τ])
  plt.plot(τ, energy)

<matplotlib.lines.Line2D

at

0x7efc79f09430>

  %%space plot
  plt.plot(τ, np.trace(rho_τ.T).real)

<matplotlib.lines.Line2D

at

0x7efc8a754ca0>

Energy Flow

  ψ_1.shape

160

500

2

Let's look at the norm.

  plt.plot(τ, (ψ_1[0].conj() * ψ_1[0]).sum(axis=1).real)

<matplotlib.lines.Line2D

at

0x7efc79f1a250>

And try to calculate the energy flow.

  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-6)
      a = L @ ψ_0.T

      return (
          eta_dot * np.array(2 * (1j * np.sum(ψ_0.conj() * a, axis=1))).flatten()
      ).real

Now we calculate the average over all trajectories.

  j = np.zeros_like(τ)
  for i in range(0, N):
      j += flow_for_traj(ψ_0[i], ψ_1[i])
  j /= N

And do the same with the alternative implementation.

  ja = np.zeros_like(τ)
  for i in range(0, N):
      ja += flow_for_traj_alt(ψ_0[i], y[i])
  ja /= N

And plot it :)

  %matplotlib inline
  plt.plot(τ, j)
  #plt.plot(τ, ja)
  plt.show()

Let's calculate the integrated energy.

  E_t = np.array([0] + [scipy.integrate.simpson(j[0:n], τ[0:n]) for n in range(1, len(τ))])
  E_t[-1]

1.9996017030561626

With this we can retrieve the energy of the interaction Hamiltonian.

  E_I = 2 - energy - E_t

  %matplotlib tk
  #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()

#+RESULTS: