Merge pull request #5 from BrieucLD/master

Add the pure dephasing at 0K example

examples/puredephasing.jl

@@ -0,0 +1,101 @@
+#=
+    Example of a Pure Dephasing Model at zero temperature and with an arbitrary temperature.  with an hard cut-off Ohmic spectral density J(ω) = 2αω when ω < ωc and 0 otherwise#
+
+    The dynamics is simulated using the T-TEDOPA method that maps the normal modes environment into a non-uniform tight-binding chain.
+
+    H = \frac{ΔE}{2} σ_z +  \frac{σ_z}{2} c_0 (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, QuadGK
+
+#----------------------------
+# Physical parameters
+#----------------------------
+
+ΔE = 0.008 # Energy of the electronic states
+
+d = 5 # number of Fock states of the chain modes
+
+N = 30 # length of the chain
+
+α = 0.01 # coupling strength
+
+s = 1 # ohmicity
+
+ωc = 0.035 # Cut-off of the spectral density J(ω)
+
+cpars = chaincoeffs_ohmic(N, α, s; ωc=ωc) # chain parameters, i.e. on-site energies ϵ_i, hopping energies t_i, and system-chain coupling c_0
+
+#-----------------------
+# Simulation parameters
+#-----------------------
+
+dt = 1.0 # time step
+
+tfinal = 300.0 # simulation time
+
+method = :TDVP1 # time-evolution method
+
+D = 6 # MPS bond dimension
+
+#---------------------------
+# MPO and initial state MPS
+#---------------------------
+
+H = puredephasingmpo(ΔE, d, N, cpars)
+
+# Initial electronic system in a superposition of 1 and 2
+ψ = zeros(2)
+ψ[1] = 1/sqrt(2)
+ψ[2] = 1/sqrt(2)
+
+
+A = productstatemps(physdims(H), state=[ψ, fill(unitcol(1,d), N)...]) # MPS representation of |ψ>|Vacuum>
+
+#---------------------------
+# Definition of observables
+#---------------------------
+
+ob1 = OneSiteObservable("sz", sz, 1)
+
+
+#-------------
+# Simulation
+#------------
+
+A, dat = runsim(dt, tfinal, A, H;
+                name = "pure dephasing model at zero temperature",
+                method = method,
+                obs = [ob1],
+                convobs = [ob1],
+                params = @LogParams(ΔE, N, d, α, s),
+                convparams = D,
+                reduceddensity=true,
+                verbose = false,
+                save = true,
+                plot = true,
+                );
+
+
+
+#----------
+# Analytical results at zero temperature 
+# (see The Theory of Open Quantum System, H.-P. Breuer & F. Petruccione 2002, Chapter 4)
+#----------
+
+Johmic(ω,s) = (2*α*ω^s)/(ωc^(s-1))
+f(ω,t) = (1 - cos(ω*t))/ω^2
+time_analytical = LinRange(0.0,tfinal,Int(tfinal))
+
+Γohmic(t) = - quadgk(x -> Johmic(x,s)*f(x,t), 0, ωc)[1]
+Decoherence_ohmic(t) = 0.5*exp(Γohmic(t))
+ListDecoherence_ohmic = [Decoherence_ohmic(t) for t in time_analytical]
+
+#-------------
+# Plots
+#------------
+
+ρ12 = sqrt.(real(dat["data/Reduced ρ"][1,2,:]).^2 .+ imag(dat["data/Reduced ρ"][1,2,:]).^2 )
+
+(plot(time_analytical, ListDecoherence_ohmic1, label="Analytics", title=L"Pure Dephasing, Ohmic $s=%$s \, ,\, T=0K$", linecolor =:black, xlabel="Time (arb. units)",ylabel="Coherence Amplitude", linewidth=4, titlefontsize=16, legend=:best, legendfont=16, xguidefontsize=16, yguidefontsize=16, tickfontsize=10))#,ylims=(0.25,0.5)))#,xticks=(0:2.5:10))))
+display(plot!(dat["data/times"],ρ12,lw=4,ls=:dash,label="Numerics"))

examples/sbm_zero_temperature.jl

@@ -1,10 +1,10 @@
-"""
+#=
     Example of a zero-temperature Spin-Boson Model with an hard cut-off Ohmic spectral density J(ω) = 2αω when ω < ωc and 0 otherwise
 
     The dynamics is simulated using the T-TEDOPA method that maps the normal modes environment into a non-uniform tight-binding chain.
 
     H = \\frac{ω_0}{2} σ_z + Δ σ_x + 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
 

src/MPSDynamics.jl

@@ -153,7 +153,7 @@ end
 
 export sz, sx, sy, numb, crea, anih, unitcol, unitrow, unitmat, spinSX, spinSY, spinSZ, SZ, SX, SY
 
-export chaincoeffs_ohmic, spinbosonmpo, methylbluempo, methylbluempo_correlated, methylbluempo_correlated_nocoupling, methylbluempo_nocoupling, ibmmpo, methylblue_S1_mpo, methylbluempo2, twobathspinmpo, xyzmpo
+export chaincoeffs_ohmic, spinbosonmpo, methylbluempo, methylbluempo_correlated, methylbluempo_correlated_nocoupling, methylbluempo_nocoupling, ibmmpo, methylblue_S1_mpo, methylbluempo2, twobathspinmpo, xyzmpo, puredephasingmpo
 
 export productstatemps, physdims, randmps, bonddims, elementmps
 
@@ -173,5 +173,7 @@ export @LogParams
 
 export MPOtoVector, MPStoVector
 
+export rhoreduced_2sites, rhoreduced_1site
+
 end
 

src/measure.jl

@@ -1019,3 +1019,57 @@ function leftcontractmps(A, O::Vector, N::Int=length(A))
     end
     return ρ
 end
+
+"""
+     rhoreduced_1site(A::Vector, site::Int=1)
+
+Caculate the reduced density matrix of the MPS A at the specified site.
+
+"""
+
+function rhoreduced_1site(A::Vector, site::Int=1)
+    N = length(A)
+    ρR = Vector{Any}(undef, N-site+1)
+    ρL = Vector{Any}(undef, site)
+    ρR[1] = ones(ComplexF64,1,1)
+    ρL[1] = ones(ComplexF64,1,1)
+    for i=N:-1:(site+1) # Build the right block, compressing the chain, from right ot left (indir=2)
+                ρR[N-i+2]= rhoAAstar(ρR[N-i+1], A[i], 2,0)
+    end
+    for i=1:(site-1)
+        ρL[i+1]= rhoAAstar(ρL[i], A[i], 1,0)
+    end
+    # Compress final virtual bondimension 
+    @tensoropt ρreduced[a,b,s,s'] := ρR[N-site+1][a0,b0] * conj(A[site][a,a0,s']) * A[site][b,b0,s]
+    @tensoropt ρreduced2[s,s'] := ρL[site][a0,b0] * ρreduced[a0,b0,s,s']
+    return ρreduced2
+end
+
+"""
+     rhoreduced_2sites(A::Vector, site::Tuple{Int, Int})
+
+Caculate the reduced density matrix of the MPS A of two neigbour sites. The resulting dimensions will be the four physical dimensions in total, 
+corresponding to the dimensions  of the two sites
+
+"""
+
+function rhoreduced_2sites(A::Vector, sites::Tuple{Int, Int})
+    N = length(A)
+    site1, site2=sites
+    ρR = Vector{Any}(undef, N-site2+1)
+    ρL = Vector{Any}(undef, site1)
+    ρR[1] = ones(ComplexF64,1,1)
+    ρL[1] = ones(ComplexF64,1,1)
+    for i=N:-1:(site2+1) # Build the right block, compressing the chain, from right ot left (indir=2)
+                ρR[N-i+2]= rhoAAstar(ρR[N-i+1], A[i], 2,0)
+    end
+    for i=1:(site1-1)
+        ρL[i+1]= rhoAAstar(ρL[i], A[i], 1,0)
+    end
+    # Compress final virtual bondimension 
+    @tensoropt ρreduced1[a,b,s,s'] := ρR[N-site2+1][a0,b0] * conj(A[site2][a,a0,s']) * A[site2][b,b0,s]
+    @tensoropt ρreduced2[a,b,s,s'] := ρL[site1][a0,b0] * conj(A[site1][a0,a,s']) * A[site1][b0,b,s]
+    @tensoropt ρreduced[s,d1,s',d2] := ρreduced2[a0,b0,d1,d2] * ρreduced1[a0,b0,s,s']
+    return ρreduced
+end
+

src/models.jl

@@ -571,3 +571,18 @@ function nearestneighbourmpo(tree_::Tree, h0, A, Ad = A')
     Ms[hn] = Ms[hn][D:D, fill(:,nc+2)...]
     return TreeNetwork(tree, Ms)
 end
+
+function puredephasingmpo(ΔE, dchain, Nchain, chainparams; tree=false)
+    u = unitmat(2)
+
+    cparams = only(chainparams[3])
+
+    Hs = (ΔE/2)*sz
+
+    M=zeros(1,3,2,2)
+    M[1, :, :, :] = up(Hs, (cparams/2)*sz, u)
+
+    chain = hbathchain(Nchain, dchain, chainparams; tree=false, reverse=false, coupletox=true)
+    return Any[M, chain...]
+end
+

src/run_1TDVP.jl

@@ -1,4 +1,4 @@
-function run_1TDVP(dt, tmax, A, H, Dmax; obs=[], timed=false, kwargs...)
+function run_1TDVP(dt, tmax, A, H, Dmax; obs=[], timed=false, reduceddensity=false, kwargs...)
     A0=deepcopy(A)
     data = Dict{String,Any}()
 
@@ -11,6 +11,10 @@ function run_1TDVP(dt, tmax, A, H, Dmax; obs=[], timed=false, kwargs...)
     for i=1:length(obs)
         push!(data, obs[i].name => reshape(exp[i], size(exp[i])..., 1))
     end
+    if reduceddensity
+        exprho = rhoreduced_1site(A0,1)
+        push!(data, "Reduced ρ" => reshape(exprho, size(exprho)..., 1))
+    end
 
     timed && (ttdvp = Vector{Float64}(undef, numsteps))
 
@@ -31,6 +35,10 @@ function run_1TDVP(dt, tmax, A, H, Dmax; obs=[], timed=false, kwargs...)
         for (i, ob) in enumerate(obs)
             data[ob.name] = cat(data[ob.name], exp[i]; dims=ndims(exp[i])+1)
         end
+        if reduceddensity
+            exprho = rhoreduced_1site(A0,1)
+            data["Reduced ρ"] = cat(data["Reduced ρ"], exprho; dims=ndims(exprho)+1)
+        end
     end
     timed && push!(data, "deltat"=>ttdvp)
     push!(data, "times" => times)