Parity

tutorials-v5/heom/heom-1e-spin-bath-model-pure-dephasing.md

@@ -5,7 +5,7 @@ jupytext:
     extension: .md
     format_name: myst
     format_version: 0.13
-    jupytext_version: 1.14.5
+    jupytext_version: 1.15.2
 kernelspec:
   display_name: Python 3 (ipykernel)
   language: python
@@ -66,7 +66,7 @@ Note that in the above, and the following, we set $\hbar = k_\mathrm{B} = 1$.
 
 ## Setup
 
-```{code-cell} ipython3
+```{code-cell}
 import contextlib
 import time
 
@@ -97,7 +97,7 @@ from qutip.solver.heom import (
 
 Let's define some helper functions for calculating correlation function expansions, plotting results and timing how long operations take:
 
-```{code-cell} ipython3
+```{code-cell}
 def cot(x):
     """ Vectorized cotangent of x. """
     return 1. / np.tan(x)
@@ -108,7 +108,7 @@ def coth(x):
     return 1. / np.tanh(x)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def plot_result_expectations(plots, axes=None):
     """ Plot the expectation values of operators as functions of time.
 
@@ -140,7 +140,7 @@ def plot_result_expectations(plots, axes=None):
     return fig
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 @contextlib.contextmanager
 def timer(label):
     """ Simple utility for timing functions:
@@ -154,7 +154,7 @@ def timer(label):
     print(f"{label}: {end - start}")
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Solver options:
 
 options = {
@@ -175,14 +175,14 @@ And let us set up the system Hamiltonian, bath and system measurement operators:
 
 Here we set $H_{sys}=0$, which means the interaction Hamiltonian and the system Hamiltonian commute, and we can compare the numerical results to a known analytical one.  We could in principle keep $\epsilon \neq 0$, but it just introduces fast system oscillations, so it is more convenient to set it to zero.
 
-```{code-cell} ipython3
+```{code-cell}
 # Defining the system Hamiltonian
 eps = 0.0  # Energy of the 2-level system.
 Del = 0.0  # Tunnelling term
 Hsys = 0.5 * eps * sigmaz() + 0.5 * Del * sigmax()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # System-bath coupling (Drude-Lorentz spectral density)
 Q = sigmaz()  # coupling operator
 
@@ -203,7 +203,7 @@ Nk = 3
 tlist = np.linspace(0, 50, 1000)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Define some operators with which we will measure the system
 # 1,1 element of density matrix - corresonding to groundstate
 P11p = basis(2, 0) * basis(2, 0).dag()
@@ -214,15 +214,15 @@ P12p = basis(2, 0) * basis(2, 1).dag()
 
 To get a non-trivial result we prepare the initial state in a superposition, and see how the bath destroys the coherence.
 
-```{code-cell} ipython3
+```{code-cell}
 # Initial state of the system.
 psi = (basis(2, 0) + basis(2, 1)).unit()
 rho0 = psi * psi.dag()
 ```
 
 ## Simulation 1: Matsubara decomposition, not using Ishizaki-Tanimura terminator
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     bath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
     HEOMMats = HEOMSolver(Hsys, bath, NC, options=options)
@@ -231,7 +231,7 @@ with timer("ODE solver time"):
     resultMats = HEOMMats.run(rho0, tlist)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Plot the results so far
 plot_result_expectations([
     (resultMats, P11p, 'b', "P11 Matsubara"),
@@ -241,7 +241,7 @@ plot_result_expectations([
 
 ## Simulation 2: Matsubara decomposition (including terminator)
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     bath = DrudeLorentzBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
     _, terminator = bath.terminator()
@@ -252,7 +252,7 @@ with timer("ODE solver time"):
     resultMatsT = HEOMMatsT.run(rho0, tlist)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Plot the results
 plot_result_expectations([
     (resultMats, P11p, 'b', "P11 Matsubara"),
@@ -266,7 +266,7 @@ plot_result_expectations([
 
 As in example 1a, we can compare to Pade and Fitting approaches.
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     bath = DrudeLorentzPadeBath(Q, lam=lam, gamma=gamma, T=T, Nk=Nk)
     HEOMPade = HEOMSolver(Hsys, bath, NC, options=options)
@@ -275,7 +275,7 @@ with timer("ODE solver time"):
     resultPade = HEOMPade.run(rho0, tlist)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Plot the results
 plot_result_expectations([
     (resultMatsT, P11p, 'b', "P11 Matsubara (+term)"),
@@ -287,7 +287,7 @@ plot_result_expectations([
 
 ## Simulation 4: Fitting approach
 
-```{code-cell} ipython3
+```{code-cell}
 def c(t, Nk):
     """ Calculates real and imag. parts of the correlation function
         using Nk Matsubara terms.
@@ -315,7 +315,7 @@ corr_ana = c(tlist_fit, lmaxmats)
 corrRana, corrIana = np.real(corr_ana), np.imag(corr_ana)
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 def wrapper_fit_func(x, N, *args):
     """ Wrapper for fitting function. """
     a, b = args[0][:N], args[0][N:2*N]
@@ -383,7 +383,7 @@ for i in range(1):
     plt.show()
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 # Set the exponential coefficients from the fit parameters
 
 ckAR = popt1[-1][0]
@@ -399,7 +399,7 @@ vkAI = -1 * popt2[-1][1]
 # vkAI = [complex(gamma)]
 ```
 
-```{code-cell} ipython3
+```{code-cell}
 with timer("RHS construction time"):
     bath = BosonicBath(Q, ckAR, vkAR, ckAI, vkAI)
     HEOMFit = HEOMSolver(Hsys, bath, NC, options=options)
@@ -410,7 +410,7 @@ with timer("ODE solver time"):
 
 ## Analytic calculations
 
-```{code-cell} ipython3
+```{code-cell}
 def pure_dephasing_evolution_analytical(tlist, wq, ck, vk):
     """
     Computes the propagating function appearing in the pure dephasing model.
@@ -487,7 +487,7 @@ def correlation_integral(t, ck, vk):
 
 For the pure dephasing analytics, we just sum up as many matsubara terms as we can:
 
-```{code-cell} ipython3
+```{code-cell}
 lmaxmats2 = 15000
 
 vk = [complex(-gamma)]
@@ -509,7 +509,7 @@ P12_ana = 0.5 * pure_dephasing_evolution_analytical(
 
 Alternatively, we can just do the integral of the propagator directly, without using the correlation functions at all
 
-```{code-cell} ipython3
+```{code-cell}
 def JDL(omega, lamc, omega_c):
     return 2. * lamc * omega * omega_c / (omega_c**2 + omega**2)
 
@@ -532,7 +532,7 @@ P12_ana2 = [
 
 ## Compare results
 
-```{code-cell} ipython3
+```{code-cell}
 plot_result_expectations([
     (resultMats, P12p, 'r', "P12 Mats"),
     (resultMatsT, P12p, 'r--', "P12 Mats + Term"),
@@ -545,7 +545,7 @@ plot_result_expectations([
 
 We can't see much difference in the plot above, so let's do a log plot instead:
 
-```{code-cell} ipython3
+```{code-cell}
 fig, axes = plt.subplots(1, 1, sharex=True, figsize=(8, 8))
 
 plot_result_expectations([
@@ -563,15 +563,15 @@ axes.legend(loc=0, fontsize=12);
 
 ## About
 
-```{code-cell} ipython3
+```{code-cell}
 qutip.about()
 ```
 
 ## Testing
 
 This section can include some tests to verify that the expected outputs are generated within the notebook. We put this section at the end of the notebook, so it's not interfering with the user experience. Please, define the tests using assert, so that the cell execution fails if a wrong output is generated.
 
-```{code-cell} ipython3
+```{code-cell}
 assert np.allclose(
     expect(P12p, resultMats.states[:15]), np.real(P12_ana)[:15],
     rtol=1e-2,