From 8b66f3833bfbc5452c2faf5903ddae22c85be7fc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Brieuc=20Le=20D=C3=A9?= Date: Tue, 5 Mar 2024 11:00:24 +0100 Subject: [PATCH] Add the pure dephasing at 0K example --- examples/puredephasing.jl | 101 +++++++++++++++++++++++++++++++ examples/sbm_zero_temperature.jl | 4 +- src/MPSDynamics.jl | 4 +- src/measure.jl | 54 +++++++++++++++++ src/models.jl | 15 +++++ src/run_1TDVP.jl | 10 ++- 6 files changed, 184 insertions(+), 4 deletions(-) create mode 100644 examples/puredephasing.jl diff --git a/examples/puredephasing.jl b/examples/puredephasing.jl new file mode 100644 index 0000000..5b029b3 --- /dev/null +++ b/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")) \ No newline at end of file diff --git a/examples/sbm_zero_temperature.jl b/examples/sbm_zero_temperature.jl index bf8bdc9..4157244 100644 --- a/examples/sbm_zero_temperature.jl +++ b/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 diff --git a/src/MPSDynamics.jl b/src/MPSDynamics.jl index e30c227..4d78756 100644 --- a/src/MPSDynamics.jl +++ b/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 diff --git a/src/measure.jl b/src/measure.jl index 4fe6aaa..7321111 100644 --- a/src/measure.jl +++ b/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 + diff --git a/src/models.jl b/src/models.jl index 89a4a6f..3191538 100644 --- a/src/models.jl +++ b/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 + diff --git a/src/run_1TDVP.jl b/src/run_1TDVP.jl index b43cfe8..64b9c91 100644 --- a/src/run_1TDVP.jl +++ b/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)