shareloqs · tfmlaX · Nov 25, 2024 · Nov 25, 2024
diff --git a/docs/src/examples/timedep.md b/docs/src/examples/timedep.md
@@ -10,10 +10,10 @@ To simulate a drive or a laser pulse, a time-dependent Hamiltonian is often need
 with ``\epsilon = \frac{2 \pi}{T_\text{Rabi}}`` the amplitude of the drive, ``T_\text{Rabi}`` the Rabi period, and ``\omega_\text{drive}`` the frequency of the drive. The drive is set up to be on resonance with the two-level system.
 
 ## The code
-First we load the `MPSdynamics.jl` package to be able to perform the simulation, the `Plots.jl` one to plot the results, and the `LaTeXStrings.jl` one to be able to use ``\LaTeX`` in the plots. The function [`MPSDynamics.disp`](@ref) is also imported.
+First we load the `MPSdynamics.jl` package to be able to perform the simulation, the `Plots.jl` one to plot the results, the `LinearAlgebra.jl` one to perform matrix diagonalization and the `LaTeXStrings.jl` one to be able to use ``\LaTeX`` in the plots. The function [`MPSDynamics.disp`](@ref) is also imported.
 
 ```julia
-using MPSDynamics, Plots, LaTeXStrings
+using MPSDynamics, Plots, LaTeXStrings, LinearAlgebra
 
 import MPSDynamics: disp
 ```
@@ -24,8 +24,6 @@ Among these, three are convergence parameters:
 *  `d` is the number of states we retain for the truncated harmonic oscillators representation of environmental modes
 * `N` is the number of chain (environmental) modes we keep. This parameters determines the maximum simulation time of the simulation: indeed excitations that arrive at the end of the chain are reflected towards the system and can lead to unphysical results
 
-The variable `Ndrive` represents the site of the MPO where the operator of the time-dependent part acts. For this example, the two-level system is at the first site of the MPS. 
-
 ```julia
 #----------------------------
 # Physical parameters
@@ -50,8 +48,6 @@ Trabi = 30.0 # Rabi period of the drive
 ϵ = 2*pi / Trabi # Intensity of the drive
 
 ωdrive = ω0 # Frequency of the drive
-
-Ndrive = 1 #Number of the site on which the drive is applied
 ```
 We set the simulation parameters and choose a time evolution method.
 As always for simulations of dynamics, the time step must be chosen wisely. The error of the TDVP methods is ``\mathcal{O}(dt^3)``.
@@ -74,6 +70,61 @@ method = :TDVP1 # time-evolution method
 
 D = [6] # MPS bond dimension
 ```
+We then set the time-dependent MPO. The variable `Ndrive` represents the site of the MPO where the operator of the time-dependent part acts. For this example, the two-level system is at the first site of the MPS. 
+
+```julia
+#---------------------------
+# MPO time-dependent
+#---------------------------
+
+timelist = collect(0:dt:tfinal)
+numsteps = length(timelist)-1
+
+#### #=
+    #Example of a Ht = σ_x ϵ sin(ωdrive t) drive on the TLS
+    #To comment if the example of (a+a^\dagger) drive is used
+
+Ndrive = 1 #Number of the site on which the drive is applied. Here Ndrive=1 for the TLS
+Ht = [ϵ*sx*sin(ωdrive*tstep) for tstep in timelist] # Time-dependent Hamiltonian term
+
+#### =#
+```
+The drive Hamiltonian can also be applied on several sites. The example of a Hamiltonian of the type :
+```math
+        \hat{H}_\text{drive}(t) = \epsilon (\hat{a}_{Ndrive} + \hat{a_{Ndrive}^\dagger}) \sin(\omega_\text{drive} t)
+```
+is illustrated in the following commented section. In order to drive one environment frequency, the operators have to be translated into chain operators. To try this type of driving, the previous section has to be commented.
+```julia
+#=
+   # Example of a Ht = (a_{Ndrive_star}+a_{Ndrive_star}^\dagger) ϵ sin(ωdrive t) drive on the (Ndrive_star)th mode. 
+   # The operators have to be transformed to drive the chain modes instead of the mode of the initial Hamiltonian
+   # This example assumes the same dimensions for all Ndrive and that applies to a chain at the right of the system.
+   # If it is not the case, the run_1TDVP loop involved when timedep=true has to be modified
+   # To comment if the example of σ_x drive is used
+
+Ndrive_star = 10
+Ndrive = collect(2:N+1)
+
+# Construct the chain Hamiltonian to write the star - chain transition matrix 
+t_chain=cpars[2][1:N-1] ; e_chain=cpars[1][1:N]
+hmat_chain = MPSDynamics.diagm(0=>e_chain, 1=>t_chain, -1=>t_chain)
+U_chain  = eigen(hmat_chain).vectors
+
+print("\n Driving of the mode of frequency : ",MPSDynamics.eigenchain(cpars, nummodes=N).values[Ndrive_star])
+
+Ht = Vector{Any}(undef,N+1)
+
+# Construct the chain operators that correspond to star operators for (a_{Ndrive_star}+a_{Ndrive_star}^\dagger
+for i in Ndrive
+    local anih_Ndrive = (U_chain[i-1,Ndrive_star]*anih(d))
+    local crea_Ndrive = (conj(U_chain[i-1,Ndrive_star])*crea(d))
+    local Ht_Ndrive = [ϵ*(anih_Ndrive+crea_Ndrive)*sin(ωdrive*tstep) for tstep in timelist]
+    Ht[i] = Ht_Ndrive
+end
+
+=#
+```
+
 Using `MPSDynamics.jl` built-in methods we define the SBM MPO and the MPS representing the initial state.
 This initial state is a product state between the system and the chain. It is constructed using a list of the 'local state' of each site of the MPS, and the dimensions of the physical legs of the MPS are set to be the same as the ones of the MPO.
 
@@ -85,11 +136,6 @@ In this part, the time-dependent terms of the MPO are stored in a list. This ena
 # MPO and initial state MPS
 #---------------------------
 
-timelist = collect(0:dt:tfinal)
-numsteps = length(timelist)-1
-
-Ht = [ϵ*sx*sin(ωdrive*tstep) for tstep in timelist] # Time-dependent Hamiltonian term
-
 H = spinbosonmpo(ω0, Δ, d, N, cpars) # MPO representation of the Hamiltonian
 
 ψ = unitcol(2,2) # Initial down-z system state 
@@ -104,7 +150,7 @@ We then choose the observables that will be stored in the data and the [`MPSDyna
 
 ob1 = OneSiteObservable("sz", sz, 1)
 
-ob2 = OneSiteObservable("chain mode occupation", numb(d), (2,N+1))
+ob2 = TwoSiteObservable("cdagc", crea(d), anih(d), collect(2:N+1), collect(2:N+1))
 
 ob3 = TwoSiteObservable("SXdisp", sx, disp(d), [1], collect(2:N+1))
 
@@ -130,13 +176,19 @@ A, dat = runsim(dt, tfinal, A, H;
                 plot = true,
                 );
 ```
-Eventually, the stored observables can be represented
+Eventually, the stored observables can be represented. For more information about the chain observables, see [Inspecting the bath by undoing the chain mapping] 
 ```julia
 
 #----------
 # Plots
 #----------
 
 plot(dat["data/times"], dat["data/sz"], label="Dmax = $(D...)", xlabel=L"t",ylabel=L"\sigma_z", title="")
+
+bath_occup = mapslices(X -> MPSDynamics.measuremodes(X, cpars[1], cpars[2]), dat["data/cdagc"], dims = [1,2])
+omeg = MPSDynamics.eigenchain(cpars, nummodes=N).values
+
+plot(omeg, bath_occup[:, :, end], lw=4, xlabel=L"\omega", ylabel=L"\langle n^b_\omega \rangle",
+title="Mode occupation in the extended bath at final time")
 ```
 
diff --git a/examples/sbm_Htimedependent.jl b/examples/sbm_Htimedependent.jl
@@ -6,7 +6,7 @@
     H = \\frac{ω_0}{2} σ_z + Δ σ_x + σ_x ϵ sin(ωdrive t) + c_0 σ_x(b_0^\\dagger + b_0) + \\sum_{i=0}^{N-1} t_i (b_{i+1}^\\dagger b_i +h.c.) + \\sum_{i=0}^{N-1} ϵ_i b_i^\\dagger b_i  
 =#
 
-using MPSDynamics, Plots, LaTeXStrings
+using MPSDynamics, Plots, LaTeXStrings, LinearAlgebra
 
 import MPSDynamics: disp
 
@@ -34,8 +34,6 @@ Trabi = 30.0 # Rabi period of the drive
 
 ωdrive = ω0 # Frequency of the drive
 
-Ndrive = 1 #Number of the site on which the drive is applied
-
 #-----------------------
 # Simulation parameters
 #-----------------------
@@ -49,14 +47,53 @@ method = :TDVP1 # time-evolution method
 D = [6] # MPS bond dimension
 
 #---------------------------
-# MPO and initial state MPS
+# MPO time-dependent
 #---------------------------
 
 timelist = collect(0:dt:tfinal)
 numsteps = length(timelist)-1
 
+#### #=
+    #Example of a Ht = σ_x ϵ sin(ωdrive t) drive on the TLS
+    #To comment if the example of (a+a^\dagger) drive is used
+
+Ndrive = 1 #Number of the site on which the drive is applied. Here Ndrive=1 for the TLS
 Ht = [ϵ*sx*sin(ωdrive*tstep) for tstep in timelist] # Time-dependent Hamiltonian term
 
+#### =#
+
+#=
+   # Example of a Ht = (a_{Ndrive_star}+a_{Ndrive_star}^\dagger) ϵ sin(ωdrive t) drive on the (Ndrive_star)th mode. 
+   # The operators have to be transformed to drive the chain modes instead of the mode of the initial Hamiltonian
+   # This example assumes the same dimensions for all Ndrive and that applies to a chain at the right of the system.
+   # If it is not the case, the run_1TDVP loop involved when timedep=true has to be modified
+   # To comment if the example of σ_x drive is used
+
+Ndrive_star = 10
+Ndrive = collect(2:N+1)
+
+# Construct the chain Hamiltonian to write the star - chain transition matrix 
+t_chain=cpars[2][1:N-1] ; e_chain=cpars[1][1:N]
+hmat_chain = MPSDynamics.diagm(0=>e_chain, 1=>t_chain, -1=>t_chain)
+U_chain  = eigen(hmat_chain).vectors
+
+print("\n Driving of the mode of frequency : ",MPSDynamics.eigenchain(cpars, nummodes=N).values[Ndrive_star])
+
+Ht = Vector{Any}(undef,N+1)
+
+# Construct the chain operators that correspond to star operators for (a_{Ndrive_star}+a_{Ndrive_star}^\dagger
+for i in Ndrive
+    local anih_Ndrive = (U_chain[i-1,Ndrive_star]*anih(d))
+    local crea_Ndrive = (conj(U_chain[i-1,Ndrive_star])*crea(d))
+    local Ht_Ndrive = [ϵ*(anih_Ndrive+crea_Ndrive)*sin(ωdrive*tstep) for tstep in timelist]
+    Ht[i] = Ht_Ndrive
+end
+
+=#
+#---------------------------
+# MPO and initial state MPS
+#---------------------------
+
 H = spinbosonmpo(ω0, Δ, d, N, cpars) # MPO representation of the Hamiltonian
 
 ψ = unitcol(2,2) # Initial down-z system state 
@@ -69,23 +106,24 @@ A = productstatemps(physdims(H), state=[ψ, fill(unitcol(1,d), N)...]) # MPS rep
 
 ob1 = OneSiteObservable("sz", sz, 1)
 
-ob2 = OneSiteObservable("chain mode occupation", numb(d), (2,N+1))
+ob2 = TwoSiteObservable("cdagc", crea(d), anih(d), collect(2:N+1), collect(2:N+1))
 
 ob3 = TwoSiteObservable("SXdisp", sx, disp(d), [1], collect(2:N+1))
 
+
 #-------------
 # Simulation
 #------------
 
 A, dat = runsim(dt, tfinal, A, H;
                 name = "Driving field on ohmic spin boson model",
                 method = method,
-                obs = [ob1],
+                obs = [ob1, ob2, ob3],
                 convobs = [ob1],
                 params = @LogParams(N, d, α, Δ, ω0, s),
                 convparams = D,
                 timedep = true, # the Hamiltonian is time dependent
-                Ndrive = Ndrive, # the first site of the MPS/MPO (i.e. the system) is concerned
+                Ndrive = Ndrive, # site(s) of the MPS/MPO that are driven
                 Htime = Ht, # list of time-dependent terms
                 verbose = false,
                 save = true,
@@ -97,3 +135,9 @@ A, dat = runsim(dt, tfinal, A, H;
 #----------
 
 plot(dat["data/times"], dat["data/sz"], label="Dmax = $(D...)", xlabel=L"t",ylabel=L"\sigma_z", title="")
+
+bath_occup = mapslices(X -> MPSDynamics.measuremodes(X, cpars[1], cpars[2]), dat["data/cdagc"], dims = [1,2])
+omeg = MPSDynamics.eigenchain(cpars, nummodes=N).values
+
+plot(omeg, bath_occup[:, :, end], lw=4, xlabel=L"\omega", ylabel=L"\langle n^b_\omega \rangle",
+title="Mode occupation in the extended bath at final time")
diff --git a/src/logging.jl b/src/logging.jl
@@ -37,11 +37,11 @@ function open_log(dt, tmax, convparams, method, machine, savedir, unid, name, pa
                 writeprintln([f,f0])
             end
             writeprintln([f,f0], "\t convparams : $convparams")
-            writeprint([f,f0], "\t options : ")
-            for (key, value) in kwargs
-                writeprint([f,f0], String(key), " = ", value, ", ")
-            end
-            writeprintln([f,f0])
+            #writeprint([f,f0], "\t options : ")
+            #for (key, value) in kwargs
+            #    writeprint([f,f0], String(key), " = ", value, ", ")
+            #end
+            #writeprintln([f,f0])
             writeprintln(f)
         end
     end

diff --git a/src/run_1TDVP.jl b/src/run_1TDVP.jl
@@ -28,7 +28,16 @@ function run_1TDVP(dt, tmax, A, H, Dmax; obs=[], timed=false, reduceddensity=fal
         if timedep
 	   Ndrive = kwargs[:Ndrive]
 	   Htime = kwargs[:Htime]
-           H0[Ndrive][1,1,:,:] = H[Ndrive][1,1,:,:] + Htime[tstep][:,:]
+           if length(Ndrive)==1
+              size(H0[Ndrive][1,1,:,:])==size(Htime[tstep][:,:]) ? nothing : throw(error("The size of Htime does not match the size of the non-interacting part of H at Ndrive"))
+              H0[Ndrive][1,1,:,:] = H[Ndrive][1,1,:,:] + Htime[tstep][:,:]
+           else
+              for i=1:length(Ndrive)
+                 site=Ndrive[i]
+                 size(H0[site][end,1,:,:])==size(Htime[site][tstep][:,:]) ? nothing : throw(error("The size of Htime does not match the size of the non-interacting part of H at Ndrive=$site"))
+                 H0[site][end,1,:,:] = H[site][end,1,:,:] + Htime[site][tstep][:,:]
+              end
+           end
         end
         if timed
             val, t, bytes, gctime, memallocs = @timed tdvp1sweep!(dt, A0, H0, F; kwargs...)

diff --git a/src/run_2TDVP.jl b/src/run_2TDVP.jl
@@ -31,8 +31,18 @@ function run_2TDVP(dt, tmax, A, H, truncerr; obs=[], Dlim=50, savebonddims=false
         if timedep
            Ndrive = kwargs[:Ndrive]
            Htime = kwargs[:Htime]
-           H0[Ndrive][1,1,:,:] = H[Ndrive][1,1,:,:] + Htime[tstep][:,:]
+           if length(Ndrive)==1
+              size(H0[Ndrive][1,1,:,:])==size(Htime[tstep][:,:]) ? nothing : throw(error("The size of Htime does not match the size of the non-interacting part of H at Ndrive"))
+              H0[Ndrive][1,1,:,:] = H[Ndrive][1,1,:,:] + Htime[tstep][:,:]
+           else
+              for i=1:length(Ndrive)
+                 site=Ndrive[i]
+                 size(H0[site][end,1,:,:])==size(Htime[site][tstep][:,:]) ? nothing : throw(error("The size of Htime does not match the size of the non-interacting part of H at Ndrive=$site"))
+                 H0[site][end,1,:,:] = H[site][end,1,:,:] + Htime[site][tstep][:,:]
+              end
+            end
         end
+
         if timed
             val, t, bytes, gctime, memallocs = @timed tdvp2sweep!(dt, A0, H0, F; truncerr=truncerr, truncdim=Dlim, kwargs...)
             println("\t","ΔT = ", t)

diff --git a/src/run_DTDVP.jl b/src/run_DTDVP.jl
@@ -34,8 +34,18 @@ function run_DTDVP(dt, tmax, A, H, prec; obs=[], effects=false, error=false, tim
         if timedep
            Ndrive = kwargs[:Ndrive]
            Htime = kwargs[:Htime]
-           H0[Ndrive][1,1,:,:] = H[Ndrive][1,1,:,:] + Htime[tstep][:,:]
+           if length(Ndrive)==1
+              size(H0[Ndrive][1,1,:,:])==size(Htime[tstep][:,:]) ? nothing : throw(error("The size of Htime does not match the size of the non-interacting part of H at Ndrive"))
+              H0[Ndrive][1,1,:,:] = H[Ndrive][1,1,:,:] + Htime[tstep][:,:]
+           else
+              for i=1:length(Ndrive)
+                 site=Ndrive[i]
+                 size(H0[site][end,1,:,:])==size(Htime[site][tstep][:,:]) ? nothing : throw(error("The size of Htime does not match the size of the non-interacting part of H at Ndrive=$site"))
+                 H0[site][end,1,:,:] = H[site][end,1,:,:] + Htime[site][tstep][:,:]
+              end
+            end
         end
+
         A0, Afull, F, info = tdvp1sweep_dynamic!(dt, A0, H0, Afull, F;
                                                  obs=obs,
                                                  prec=prec,