mirror of
https://github.com/vale981/master-thesis
synced 2025-03-14 15:16:39 -04:00
548 lines
15 KiB
Org Mode
548 lines
15 KiB
Org Mode
#+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:
|
||
| <matplotlib.lines.Line2D | at | 0x7f6798323760> |
|
||
| <matplotlib.lines.Line2D | at | 0x7f6799a6ea60> |
|
||
| <matplotlib.lines.Line2D | at | 0x7f6799a6eca0> |
|
||
[[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)
|
||
<ipython-input-777-b4dd44af7ffd> in <module>
|
||
----> 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)
|
||
<ipython-input-781-12d991d6efe0> in <module>
|
||
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 <class '__main__.GaussBCF'>: 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:
|
||
| <matplotlib.lines.Line2D | at | 0x7f67981b3e80> |
|
||
[[file:./.ob-jupyter/892241bc3127d1dc4581b7f31e8b392425c9716c.png]]
|
||
:END:
|
||
|
||
#+begin_src jupyter-python
|
||
%%space plot
|
||
plt.plot(τ, np.trace(rho_τ.T).real)
|
||
#+end_src
|
||
|
||
#+RESULTS:
|
||
:RESULTS:
|
||
| <matplotlib.lines.Line2D | at | 0x7f6799c2caf0> |
|
||
[[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:
|
||
| <matplotlib.lines.Line2D | at | 0x7f679a203700> |
|
||
[[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)
|
||
<ipython-input-787-b9081128ed40> in <module>
|
||
1 j = np.zeros_like(τ)
|
||
2 for i in range(0, N):
|
||
----> 3 j += flow_for_traj(ψ_0[i], ψ_1[i])
|
||
4 j /= N
|
||
|
||
<ipython-input-786-c57f86a6b31b> 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)
|
||
<ipython-input-788-5dc6ccf09941> in <module>
|
||
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
|
||
|
||
<ipython-input-786-c57f86a6b31b> 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:
|
||
| <matplotlib.lines.Line2D | at | 0x7f6798fb5eb0> |
|
||
| <matplotlib.lines.Line2D | at | 0x7f6798f8cb50> |
|
||
| <matplotlib.lines.Line2D | at | 0x7f6798f75d60> |
|
||
: Text(0.5, 0, 'τ')
|
||
: <matplotlib.legend.Legend at 0x7f6798f8ceb0>
|
||
[[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)
|
||
| <matplotlib.lines.Line2D | at | 0x7f678d445f40> |
|
||
[[file:./.ob-jupyter/aaa6bdcc61452758e092eccc2d62ef8cdfe5f258.png]]
|
||
:END:
|