diff --git a/python/energy_flow_proper/08_dynamic_one_bath/coupling_and_system_cycle.org b/python/energy_flow_proper/08_dynamic_one_bath/coupling_and_system_cycle.org new file mode 100644 index 0000000..f3a35d6 --- /dev/null +++ b/python/energy_flow_proper/08_dynamic_one_bath/coupling_and_system_cycle.org @@ -0,0 +1,380 @@ +#+PROPERTY: header-args :session 08_coup_sys_cycle :kernel python :pandoc no :async yes :tangle no +Now let's try to actually extract energy. + +* Boilerplate +#+begin_src jupyter-python :results none + import figsaver as fs + from hiro_models.one_qubit_model import QubitModel, StocProcTolerances + import hiro_models.model_auxiliary as aux + import numpy as np + import qutip as qt +#+end_src + +Init ray and silence stocproc. +#+begin_src jupyter-python + import ray + ray.shutdown() + #ray.init(address='auto') + ray.init() +#+end_src + +#+RESULTS: +: RayContext(dashboard_url='', python_version='3.9.10', ray_version='1.12.0', ray_commit='f18fc31c7562990955556899090f8e8656b48d2d', address_info={'node_ip_address': '141.30.17.221', 'raylet_ip_address': '141.30.17.221', 'redis_address': None, 'object_store_address': '/tmp/ray/session_2022-04-28_16-22-15_573871_76439/sockets/plasma_store', 'raylet_socket_name': '/tmp/ray/session_2022-04-28_16-22-15_573871_76439/sockets/raylet', 'webui_url': '', 'session_dir': '/tmp/ray/session_2022-04-28_16-22-15_573871_76439', 'metrics_export_port': 64264, 'gcs_address': '141.30.17.221:57854', 'address': '141.30.17.221:57854', 'node_id': '3d8ae351108358a0a4635b495c616f4e81a6fc9404e1d11c82389af6'}) + +#+begin_src jupyter-python :results none + from hops.util.logging_setup import logging_setup + import logging + logging_setup(logging.INFO) +#+end_src + +* Cycle +#+begin_src jupyter-python :results none + from hops.util.dynamic_matrix import SmoothStep, Periodic, ConstantMatrix, ScaleTime, Shift, Harmonic +#+end_src + +Now let's design an asymetric cycle. +#+begin_src jupyter-python :results none + L_op = (1/2 * qt.sigmam() + 0 * qt.sigmax()).full() + L_pos = SmoothStep(L_op, 0, .4, 3) - SmoothStep(L_op, .6, .8, 3) + L = ScaleTime(Periodic(L_pos, 1), 1) +#+end_src + +Let's plot such a nice smooth-step. +#+begin_src jupyter-python + ts = np.linspace(0, 5, 10000) + plt.plot(ts, L.operator_norm(ts)) +#+end_src + +#+RESULTS: +:RESULTS: +| | +[[file:./.ob-jupyter/caedb7196fc09f0ed309141321acb49502ee3df5.svg]] +:END: + +Now we turn the system around after each fill-cycle. +#+begin_src jupyter-python :results none + H_op = 1 * (1/2 * (qt.sigmaz() + qt.identity(2))).full() + H = ScaleTime((ConstantMatrix(H_op) - .5 * Harmonic(qt.sigmay().full(), 1, 0)), np.pi) +#+end_src + +#+begin_src jupyter-python + plt.plot(ts, L(ts)[:,1,0].real) + plt.plot(ts, H.operator_norm(ts)) +#+end_src + +#+RESULTS: +:RESULTS: +| | +[[file:./.ob-jupyter/3760d5d9641934a8fdf7a9c97bac3d154e935f47.svg]] +:END: + +* Model: Very Non-Markovian +#+begin_src jupyter-python + model = QubitModel( + δ=.5, + ω_c=.1, + t=np.linspace(0, 20, 1_000), + ψ_0=qt.basis([2], [1]), + description=f"Testing the time dependent coupling with smooth step.", + k_max=7, + bcf_terms=6, + truncation_scheme="simplex", + driving_process_tolerance=StocProcTolerances(1e-3, 1e-3), + thermal_process_tolerance=StocProcTolerances(1e-3, 1e-3), + T = 4, + L = L, + H = H, + ) +#+end_src + +#+RESULTS: + +#+begin_src jupyter-python + plt.plot(model.t, model.spectral_density(model.t)) + plt.axvline(model.ω_c) + # plt.axhline(1/(np.exp(model.ω_c / model.T) - 1) * model.ω_c) + # plt.axhline(1/(np.exp(1e-3 / model.T) - 1) * 1e-3) + #plt.axhline(model.T ** 2 * np.pi **2 /6 * model.ω_c) +#+end_src + +#+RESULTS: +:RESULTS: +: +[[file:./.ob-jupyter/93c7bcaec0d9c6529fbf419915378372b7fb7f31.svg]] +:END: + + +#+begin_src jupyter-python + aux.integrate(model, 10000) +#+end_src + +#+RESULTS: +#+begin_example + /home/hiro/src/hops/hops/util/bcf.py:280: UserWarning: this implementation uses mpmath to evaluate the zeta_function! for a better performance consider the 'OhmEnv' package + warnings.warn( + n:32 d:0.00014166421319021928 tol:0.001 + done! + [INFO hops.core.integration 76439] Choosing the nonlinear integrator. + [INFO hops.core.integration 76439] Using 4 integrators. + [INFO hops.core.integration 76439] Some 0 trajectories have to be integrated. + [INFO hops.core.integration 76439] Using 1716 hierarchy states. + 0it [00:00, ?it/s] +#+end_example + + +#+begin_src jupyter-python + import scipy + + with aux.get_data(model) as data: + fs.plot_with_σ(model.t, model.system_energy(data)) + plt.plot( + model.t, [-np.trace(ρ @ scipy.linalg.logm(ρ)) for ρ in data.rho_t_accum.mean] + ) + plt.plot(model.t, model.H(model.t)[:, 0, 0].real) + plt.plot(model.t, model.L(model.t)[:, 0, 1].real) +#+end_src + +#+RESULTS: +:RESULTS: +: /nix/store/6md8j2v3a9xhxllmyha3gn1qwrfllm3n-python3-3.9.10-env/lib/python3.9/site-packages/scipy/linalg/_matfuncs_inv_ssq.py:827: LogmExactlySingularWarning: The logm input matrix is exactly singular. +: warnings.warn(exact_singularity_msg, LogmExactlySingularWarning) +: /nix/store/6md8j2v3a9xhxllmyha3gn1qwrfllm3n-python3-3.9.10-env/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part +: return np.asarray(x, float) +[[file:./.ob-jupyter/8142cb1061478cc43e470575d8d7382628cbca4b.svg]] +:END: + + +#+begin_src jupyter-python + fig, ax = fs.plot_energy_overview(model) + plt.plot(model.t, model.H.operator_norm(model.t), label="$H(t)$") + plt.plot(model.t, model.L.operator_norm(model.t), label="$L(t)$") + with aux.get_data(model) as data: + plt.plot( + model.t, [-np.trace(ρ @ scipy.linalg.logm(ρ)) for ρ in data.rho_t_accum.mean], + label="$S[ρ_{s}]$" + ) + ax.legend() +#+end_src + +#+RESULTS: +:RESULTS: +#+begin_example + /home/hiro/src/hops/hops/util/bcf.py:280: UserWarning: this implementation uses mpmath to evaluate the zeta_function! for a better performance consider the 'OhmEnv' package + warnings.warn( + n:32 d:0.00014166421319021928 tol:0.001 + done! + n:32 d:0.00014166421319021928 tol:0.001 + done! + n:32 d:0.00014166421319021928 tol:0.001 + done! + /nix/store/6md8j2v3a9xhxllmyha3gn1qwrfllm3n-python3-3.9.10-env/lib/python3.9/site-packages/scipy/linalg/_matfuncs_inv_ssq.py:827: LogmExactlySingularWarning: The logm input matrix is exactly singular. + warnings.warn(exact_singularity_msg, LogmExactlySingularWarning) + /nix/store/6md8j2v3a9xhxllmyha3gn1qwrfllm3n-python3-3.9.10-env/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part + return np.asarray(x, float) +#+end_example +: +[[file:./.ob-jupyter/dec6c0d24496e2029b219adca1c4496a70fcbacc.svg]] +:END: + + + +#+begin_src jupyter-python :results none + with aux.get_data(model) as data: + ρ = data.rho_t_accum.mean[:] + σ_ρ = data.rho_t_accum.ensemble_std[:] + xs = np.einsum("tij,ji->t", ρ, qt.sigmax().full()) + ys = np.einsum("tij,ji->t", ρ, qt.sigmay().full()) + zs = np.einsum("tij,ji->t", ρ, qt.sigmaz().full()) +#+end_src + + +#+begin_src jupyter-python + b = qt.Bloch() + b.add_points([xs, ys, zs]) + b.view = [20, 20] + b.point_size = [0.1] + b.sphere_alpha = 0.1 + b.size = [20, 20] + b.render() + b.fig +#+end_src + +#+RESULTS: +[[file:./.ob-jupyter/ee1a474570fad49dfefa176ff8b1b3262354d107.svg]] + + +#+begin_src jupyter-python + #plt.plot(model.t, zs) + #plt.plot(model.t, xs) + plt.plot(model.t, np.abs(ρ[:, 0,1])) + plt.plot(model.t, σ_ρ[:, 0,1]) +#+end_src + +#+RESULTS: +:RESULTS: +| | +[[file:./.ob-jupyter/155688dabfd7fbccbdab22a83c9cd9f1aa34f5d1.svg]] +:END: + +* Model: Maximum Flow +#+begin_src jupyter-python + model = QubitModel( + δ=.5, + ω_c=2, + t=np.linspace(0, 5, 1_000), + ψ_0=qt.basis([2], [1]), + description=f"Testing the time dependent coupling with smooth step.", + k_max=7, + bcf_terms=6, + truncation_scheme="simplex", + driving_process_tolerance=StocProcTolerances(1e-3, 1e-3), + thermal_process_tolerance=StocProcTolerances(1e-3, 1e-3), + T = 4, + L = L, + H = ConstantMatrix(H_op), + ) +#+end_src + +#+RESULTS: + + +#+begin_src jupyter-python + aux.integrate(model, 100) +#+end_src + +#+RESULTS: +#+begin_example + /home/hiro/src/hops/hops/util/bcf.py:280: UserWarning: this implementation uses mpmath to evaluate the zeta_function! for a better performance consider the 'OhmEnv' package + warnings.warn( + n:32 d:0.33446500514383937 tol:0.001 + n:64 d:0.02262163635902933 tol:0.001 + n:128 d:2.2746881114193725e-05 tol:0.001 + done! + [INFO hops.core.integration 76439] Choosing the nonlinear integrator. + [INFO hops.core.integration 76439] Using 4 integrators. + [INFO hops.core.integration 76439] Some 100 trajectories have to be integrated. + [INFO hops.core.integration 76439] Using 1716 hierarchy states. +100% 100/100 [01:21<00:00, 1.23it/s] +#+end_example + + +#+begin_src jupyter-python + import scipy + + with aux.get_data(model) as data: + fs.plot_with_σ(model.t, model.system_energy(data)) + plt.plot( + model.t, [-np.trace(ρ @ scipy.linalg.logm(ρ)) for ρ in data.rho_t_accum.mean] + ) + plt.plot(model.t, model.H(model.t)[:, 0, 0].real) + plt.plot(model.t, model.L(model.t)[:, 0, 1].real) +#+end_src + +#+RESULTS: +:RESULTS: +: Loading: 100% 5/5 [00:00<00:00, 82.94it/s] +: Processing: 100% 5/5 [00:00<00:00, 761.02it/s] +: [INFO root 76439] Writing cache to: results/3694939998333a5d25b932d5b63ea70f47cb547297a909f142e8b51d8a096d9b_573fc383266c5ec57f5ac393c4145a64be86ba7915c0cd8da8bfc887bdb3cbc0_op_exp_task_100_None_1ff6d34b7de893995dfd0ae8efe2ed770a6773ef52b6e717b4b4ee9b9e5d285d.npy +: /nix/store/6md8j2v3a9xhxllmyha3gn1qwrfllm3n-python3-3.9.10-env/lib/python3.9/site-packages/scipy/linalg/_matfuncs_inv_ssq.py:827: LogmExactlySingularWarning: The logm input matrix is exactly singular. +: warnings.warn(exact_singularity_msg, LogmExactlySingularWarning) +: /nix/store/6md8j2v3a9xhxllmyha3gn1qwrfllm3n-python3-3.9.10-env/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part +: return np.asarray(x, float) +[[file:./.ob-jupyter/84f50c19b292b6f8968709211eb9f1e8e8b83a91.svg]] +:END: + + +#+begin_src jupyter-python + fig, ax = fs.plot_energy_overview(model) + plt.plot(model.t, model.H.operator_norm(model.t), label="$H(t)$") + plt.plot(model.t, model.L.operator_norm(model.t), label="$L(t)$") + with aux.get_data(model) as data: + plt.plot( + model.t, [-np.trace(ρ @ scipy.linalg.logm(ρ)) for ρ in data.rho_t_accum.mean], + label="$S[ρ_{s}]$" + ) + ax.legend() +#+end_src + +#+RESULTS: +:RESULTS: +#+begin_example + /home/hiro/src/hops/hops/util/bcf.py:280: UserWarning: this implementation uses mpmath to evaluate the zeta_function! for a better performance consider the 'OhmEnv' package + warnings.warn( + n:32 d:0.33446500514383937 tol:0.001 + n:64 d:0.02262163635902933 tol:0.001 + n:128 d:2.2746881114193725e-05 tol:0.001 + done! + Loading: 100% 5/5 [00:00<00:00, 19.33it/s] + Processing: 100% 5/5 [00:00<00:00, 22.99it/s] + [INFO root 76439] Writing cache to: results/flow_573fc383266c5ec57f5ac393c4145a64be86ba7915c0cd8da8bfc887bdb3cbc0_flow_worker_100_None_0c0cafd9366880235df0a72a856be2f076e35acf3353bfe08d1d3ec198a02f32.npy + n:32 d:0.33446500514383937 tol:0.001 + n:64 d:0.02262163635902933 tol:0.001 + n:128 d:2.2746881114193725e-05 tol:0.001 + done! + Loading: 100% 5/5 [00:00<00:00, 30.33it/s] + Processing: 100% 5/5 [00:00<00:00, 19.61it/s] + [INFO root 76439] Writing cache to: results/interaction_573fc383266c5ec57f5ac393c4145a64be86ba7915c0cd8da8bfc887bdb3cbc0_interaction_energy_task_100_None_a1ebae168ba6027d405a66e94ca963c84f8aedb16621da5954b5631bcd1636c4.npy + n:32 d:0.33446500514383937 tol:0.001 + n:64 d:0.02262163635902933 tol:0.001 + n:128 d:2.2746881114193725e-05 tol:0.001 + done! + /nix/store/6md8j2v3a9xhxllmyha3gn1qwrfllm3n-python3-3.9.10-env/lib/python3.9/site-packages/scipy/linalg/_matfuncs_inv_ssq.py:827: LogmExactlySingularWarning: The logm input matrix is exactly singular. + warnings.warn(exact_singularity_msg, LogmExactlySingularWarning) + /nix/store/6md8j2v3a9xhxllmyha3gn1qwrfllm3n-python3-3.9.10-env/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part + return np.asarray(x, float) +#+end_example +: +[[file:./.ob-jupyter/494a441d03c95f40fbd6b286f580057e33d5f49c.svg]] +:END: + + + +#+begin_src jupyter-python :results none + with aux.get_data(model) as data: + ρ = data.rho_t_accum.mean[:] + σ_ρ = data.rho_t_accum.ensemble_std[:] + xs = np.einsum("tij,ji->t", ρ, qt.sigmax().full()) + ys = np.einsum("tij,ji->t", ρ, qt.sigmay().full()) + zs = np.einsum("tij,ji->t", ρ, qt.sigmaz().full()) +#+end_src + + +#+begin_src jupyter-python + b = qt.Bloch() + b.add_points([xs, ys, zs]) + b.view = [20, 20] + b.point_size = [0.1] + b.sphere_alpha = 0.1 + b.size = [20, 20] + b.render() + b.fig +#+end_src + +#+RESULTS: +[[file:./.ob-jupyter/3f303eb2cb2fe0ca8714be20de2b8856aa6e6951.svg]] + + +#+begin_src jupyter-python + #plt.plot(model.t, zs) + #plt.plot(model.t, xs) + plt.plot(model.t, np.abs(ρ[:, 0,1])) + plt.plot(model.t, σ_ρ[:, 0,1]) +#+end_src + +#+RESULTS: +:RESULTS: +| | +[[file:./.ob-jupyter/00b87c2e8644894417234a73c220d39a811a5aa5.svg]] +:END: + + +* Findings +- flipping H does not work -> Thermalization does not create inversion +- one can't extract energy from the single bath because it is in a + passive state -> WRONG +- Choosing $T>2*\omega_c$ and $\omega_c\approx ||H||$ and $L=\sigma_{+}$ yields the deal +- after certain time span -> non-passivity exausted -> when system is "full" +- adding some $\sigma_x$ reduces the total intial energy loss and hampers energy transfer + - turn-around comes later because energy transfer is slower +- maybe zero mean entropy production helps? +- maybe the "forgetful" nature of the bath helps