init ergo stuff

python/energy_flow_proper/ergo_stuff/.envrc

@@ -0,0 +1,2 @@
+use_flake
+eval "$shellHook"

python/energy_flow_proper/ergo_stuff/ergotropy_bath_qubit.org

@@ -0,0 +1,277 @@
+#+PROPERTY: header-args :session 09_ergo :kernel python :pandoc no :async yes :tangle no
+
+#+begin_src jupyter-python :results none
+  import figsaver as fs
+  import numpy as np
+#+end_src
+
+* Ergotropy with Two HO
+Let's generate our initial states.
+#+begin_src jupyter-python
+  import itertools, pprint, functools
+
+  N = 100
+  ω_1 = 1
+  ω_2 = 3
+
+
+  def energy(state, ω_1, ω_2):
+      n1, n2, _ = state
+      return ω_1 * n1 + ω_2 * n2
+
+
+  def weight(state, ω_1, ω_2):
+      return np.exp(-energy(state, ω_1, ω_2)) * state[-1]
+
+
+  @functools.cache
+  def gen_states(N):
+      return [(n1, n2, 0) for n1, n2 in itertools.product(range(N), range(N))] + [
+          (n1, n2, 1) for n1, n2 in itertools.product(range(N), range(N))
+      ]
+
+
+  def order_by_energy(states, ω_1, ω_2):
+      return sorted(states, key=lambda state: energy(state, ω_1, ω_2))
+
+
+  def order_by_weight(states, ω_1, ω_2):
+      return sorted(states, key=lambda state: -weight(state, ω_1, ω_2))
+
+
+  def initial_energy(states, ω_1, ω_2):
+      return sum(energy(state, ω_1, ω_2) * weight(state, ω_1, ω_2) for state in states) / sum(weight(state, ω_1, ω_2) for state in states)
+
+
+  def reorder_states(states, ω_1, ω_2):
+      weight_ordered_states = order_by_weight(states, ω_1, ω_2)
+      energy_ordered_states = order_by_energy(states, ω_1, ω_2)
+
+      return (
+          (
+              energy_ordered_state,
+              (
+                  weight(weight_ordered_state, ω_1, ω_2),
+                  energy(energy_ordered_state, ω_1, ω_2),
+              ),
+          )
+          for energy_ordered_state, weight_ordered_state in zip(energy_ordered_states, weight_ordered_states)
+      )
+
+
+  def calculate_ergotropy(N, ω_1, ω_2):
+      states = list(gen_states(N))
+      reordered_states = list(reorder_states(states, ω_1, ω_2))
+      ergo = initial_energy(states, ω_1, ω_2) - sum(
+          state[1][0] * state[1][1] for state in reordered_states
+      ) / sum(
+          state[1][0] for state in reordered_states
+      )
+
+      return ergo
+
+  def reduced_state(N, ω_1, ω_2):
+      states = list(gen_states(N))
+      reordered_states = reorder_states(states, ω_1, ω_2)
+
+      weights = [0, 0]
+
+      for state in reordered_states:
+          weights[state[0][-1]] += state[1][0]
+
+      total = sum(weights)
+      return [w / total for w in weights]
+
+  def ergo_one_ho(ω):
+      return (np.exp(-ω) - np.exp(-2 * ω)) / (1 - np.exp(-ω) - np.exp(-2 * ω) + np.exp(-3 * ω)) * ω
+#+end_src
+
+#+RESULTS:
+
+
+Now we fill up our degeneracies.
+#+begin_src jupyter-python
+  ω_2s = np.linspace(1/2, 10, 100)
+  ergos = [calculate_ergotropy(100, 1, two) for two in ω_2s]
+  plt.plot(ω_2s, ergos)
+  plt.plot(
+      ω_2s,
+      ergo_one_ho(ω_2s),
+      label="one HO",
+  )
+
+  plt.plot(
+      ω_2s,
+      (ergo_one_ho(ω_2s) + ergo_one_ho(1))/2,
+      label="mean two single-HO",
+  )
+
+  plt.axhline(ergo_one_ho(1), label="Ergo of one HO with ω=1")
+  plt.legend()
+  plt.ylabel("Ergotropy")
+  plt.xlabel("ω")
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+# [goto error]
+: ---------------------------------------------------------------------------
+: NameError                                 Traceback (most recent call last)
+: Input In [3], in <cell line: 1>()
+: ----> 1 ω_2s = np.linspace(1/2, 10, 100)
+:       2 ergos = [calculate_ergotropy(100, 1, two) for two in ω_2s]
+:       3 plt.plot(ω_2s, ergos)
+:
+: NameError: name 'np' is not defined
+:END:
+
+For ω_2->∞ we get the ergotropy of the lowest HO. Interestingly -> not ergo of TWO ho when same Freq
+
+* Ergotropy with N HO
+Let's generate our initial states.
+#+begin_src jupyter-python
+  import itertools, pprint, functools
+
+
+  def energy(state, ω):
+      n, _ = state
+      return sum(ω * n)
+
+
+  def weight(state, ω):
+      return np.exp(-energy(state, ω)) * (1-state[-1])
+
+
+  def gen_states(ω, max_dim, min_contrib):
+      N = len(ω)
+      return list((
+          (n, i)
+          for n in itertools.product(*((range(max_dim),) * N))
+          if (energy((n, 0), ω) * weight((n, 0), ω)) >= min_contrib
+          for i in range(2)
+     ))
+
+
+  def order_by_energy(states, ω):
+      return sorted(states, key=lambda state: energy(state, ω))
+
+
+  def order_by_weight(states, ω):
+      return sorted(states, key=lambda state: -weight(state, ω))
+
+
+  def initial_energy(states, ω):
+      return sum(
+          energy(state, ω) * weight(state, ω) for state in states
+      ) / sum(weight(state, ω) for state in states)
+
+
+  def reorder_states(states, ω):
+      weight_ordered_states = order_by_weight(states, ω)
+      energy_ordered_states = order_by_energy(states, ω)
+
+      return (
+          (
+              energy_ordered_state,
+              (
+                  weight(weight_ordered_state, ω),
+                  energy(energy_ordered_state, ω),
+              ),
+          )
+          for energy_ordered_state, weight_ordered_state in zip(
+              energy_ordered_states, weight_ordered_states
+          )
+      )
+
+
+  def calculate_ergotropy(states, ω):
+      reordered_states = list(reorder_states(states, ω))
+
+      ergo = initial_energy(states, ω) - sum(
+          state[1][0] * state[1][1] for state in reordered_states
+      ) / sum(state[1][0] for state in reordered_states)
+
+      return ergo
+
+
+  def reduced_state(states, ω):
+      reordered_states = reorder_states(states, ω)
+
+      weights = [0, 0]
+
+      for state in reordered_states:
+          weights[state[0][-1]] += state[1][0]
+
+      total = sum(weights)
+      return [w / total for w in weights]
+
+#+end_src
+
+#+RESULTS:
+
+
+#+begin_src jupyter-python
+  def ergo_spread(ω_min, Δ_max, N=2, dim=3, min_weight=1e-5, n=100):
+      Δs = np.linspace(0.0, Δ_max, n)
+      ω_s = [ω_min + np.linspace(0, 1, N) * Δ for Δ in Δs]
+      ergos = [calculate_ergotropy(gen_states(ω, dim, 1e-5), ω) for ω in ω_s]
+      return Δs, ω_s, ergos
+
+
+  @fs.wrap_plot
+  def plot_ergo_spread(*args, ax=None, **kwargs):
+      Δs, ω_s, ergos = ergo_spread(*args, **kwargs)
+      ax.plot(Δs, ergos, label=fr"exact, $ω_m={ω_c}, N={len(ω_s[0])}$")
+
+      # ax.plot(
+      #     Δs,
+      #     [np.sum([ergo_one_ho(ω_i) for ω_i in ω]) for ω in ω_s],
+      #     label="sum of single-HOs",
+      # )
+
+      # ax.axhline(ergo_one_ho(np.min(ω_s)), label=f"Ergo of one HO with ω={np.min(ω_s)}")
+      ax.legend()
+      ax.set_ylabel("Ergotropy")
+      ax.set_xlabel(rf"$Δ$")
+
+      return Δs, ω_s, ergos
+#+end_src
+
+#+RESULTS:
+
+
+#+begin_src jupyter-python
+  fig, ax = plt.subplots()
+  for N in range(10, 16):
+      plot_ergo_spread(1, 3, N, 2, n=5, ax=ax)
+#+end_src
+
+#+RESULTS:
+[[file:./.ob-jupyter/8244f1653cfe1f498f95127ea7668d2c75976aff.svg]]
+
+Limit works still. But depends on the absolute difference of the first two i believe.
+Interesingly, the ergotropy drops when we reduce the spacing. There is a nontrivial maximum.
+The effect of big spacing is explained by the suppression of excitations in the high freq. oscis.
+
+There appears to be an upper bound.
+
+** TODO Verify That N-Times the same osci does not increase Ergo asymptotically (for low T)
+#+begin_src jupyter-python :results none
+  Ns = np.array(range(1,20))
+  ω_f = 2
+  ergos = [calculate_ergotropy(gen_states(np.repeat(ω_f, N), 2, 1e-8), np.repeat(ω_f, N)) for N in Ns]
+#+end_src
+
+#+begin_src jupyter-python
+  plt.plot(Ns, ergos)
+#+end_src
+
+#+RESULTS:
+:RESULTS:
+| <matplotlib.lines.Line2D | at | 0x7ff95eee8760> |
+[[file:./.ob-jupyter/77d8c3967791e403c16c04f63279e6a852edbd0a.svg]]
+:END:
+
+Hmm. for Qubit Bath it works. That should actually also work analytically.
+
+*** TODO Verify that Entropy of requced Qubit State approaches Maximum

python/energy_flow_proper/ergo_stuff/flake.nix

@@ -0,0 +1,27 @@
+{
+  description = "Expoloring the ergotropy of (finite) HO baths.";
+
+  inputs = {
+    nixpkgs.url = "nixpkgs/nixos-unstable";
+    utils.url = "github:vale981/hiro-flake-utils";
+  };
+
+  outputs = { self, utils, nixpkgs, ... }:
+    (utils.lib.poetry2nixWrapper nixpkgs {
+      name = "ergo_stuff";
+      shellPackages = pkgs: with pkgs; [ pyright python39Packages.jupyter sshfs ];
+      python = pkgs: pkgs.python39Full;
+      shellOverride = (oldAttrs: {
+        shellHook = ''
+                    # export PYTHONPATH=/home/hiro/src/two_qubit_model/:$PYTHONPATH
+                    # export PYTHONPATH=/home/hiro/src/hops/:$PYTHONPATH
+                    # export PYTHONPATH=/home/hiro/src/hopsflow/:$PYTHONPATH
+                    # export PYTHONPATH=/home/hiro/src/stocproc/:$PYTHONPATH
+                    '';
+      });
+      noPackage = true;
+      poetryArgs = {
+        projectDir = ./.;
+      };
+    });
+}

python/energy_flow_proper/ergo_stuff/pyproject.toml

@@ -0,0 +1,19 @@
+[tool.poetry]
+name = "ergo_stuff"
+version = "1.0.0"
+description = "Expoloring the ergotropy of (finite) HO baths."
+authors = ["Valentin Boettcher <hiro@protagon.space>"]
+license = "GPLv3"
+
+[tool.poetry.dependencies]
+python = ">=3.9,<3.11"
+matplotlib = "^3.5.0"
+jupyter = "^1.0.0"
+
+[tool.poetry.dev-dependencies]
+black = "^21.12b0"
+click = "==8.0.4"
+
+[build-system]
+requires = ["poetry-core>=1.0.0"]
+build-backend = "poetry.core.masonry.api"