Add puredephasing temperature example

ChainOhmT/chainOhmT.jl

@@ -1,5 +1,6 @@
 module Chaincoeffs
 using CSV
+using HDF5
 using Tables
 include("quadohmT.jl")
 include("mcdis2.jl")
@@ -10,10 +11,12 @@ Code translated in Julia by Thibaut Lacroix in 10/2020 from a previous Matlab co
 
 export mc, mp, iq, idelta, irout, AB, a, wc, beta, DM, uv
 ## Spectral density parameters
-a = 0.03
-wc = 1
-beta = 1
-xc = wc
+a = 0.01 # Coupling parameter
+wc = 0.035 # Cutoff frequency. Often set up as 1 for arbitrary units
+s = 1 #Ohmicity parameter
+beta = 100 # 1/(kB T)
+
+xc = wc # Limit of definition segments
 
 ## discretisation parameters
 mc=4 # the number of component intervals
@@ -22,7 +25,7 @@ iq=1 # a parameter to be set equal to 1, if the user provides his or her own qua
 idelta=2 # a parameter whose default value is 1, but is preferably set equal to 2, if iq=1 and the user provides Gauss-type quadrature routines
 irout=2 # choice between the Stieltjes (irout = 1) and the Lanczos procedure (irout != 1)
 AB =[[-Inf -xc];[-xc 0];[0 xc];[xc Inf]] # component intervals
-N=1000 #Number of bath modes
+N=30 #Number of bath modes
 Mmax=5000
 eps0=1e7*eps(Float64)
 
@@ -41,8 +44,26 @@ for i = 1:mc
     xw = quadfinT(Mcap,i,uv)
     global eta += sum(xw[:,2])
 end
-jacerg[N,2] = sqrt(eta/pi)
+jacerg[N,2] = sqrt(eta)
+
+# Write a CSV file
+curdir = @__DIR__
+
 header = Vector([Symbol("α_n"), Symbol("β_n and η")])
-CSV.write("chaincoeffs_ohmic_a$(a)wc$(wc)xc$(xc)beta$(beta).csv", Tables.table(jacerg, header=header))
+CSV.write("$curdir/ohmicT/chaincoeffs_ohmic_a$(a)wc$(wc)xc$(xc)beta$(beta).csv", Tables.table(jacerg, header=header))
+
+Nstr=string(N)
+astr=string(a)
+sstr=string(s)
+bstr=string(beta)
+# the "path" to the data inside of the h5 file is beta -> alpha -> s -> data (e, t or c)
+
+# Write onsite energies
+h5write("$curdir/ohmicT/chaincoeffs.h5", string("/",Nstr,"/",astr,"/",sstr,"/",bstr,"/e"), jacerg[1:N,1])
+# Write hopping energies
+h5write("$curdir/ohmicT/chaincoeffs.h5", string("/",Nstr,"/",astr,"/",sstr,"/",bstr,"/t"), jacerg[1:N-1,2])
+# Write coupling coefficient
+h5write("$curdir/ohmicT/chaincoeffs.h5", string("/",Nstr,"/",astr,"/",sstr,"/",bstr,"/c"), jacerg[N,2])
 
 end
+

ChainOhmT/ohmicT/chaincoeffs_ohmic_a0.01wc0.035xc0.035beta100.csv

@@ -0,0 +1,31 @@
+α_n,β_n and η
+0.015637470905482898,0.015392189234865787
+0.0010976697115897855,0.018118271985527956
+-0.000791513873033027,0.017940785434537584
+-0.0003291962788118087,0.01769964595571843
+-0.00012347015587491842,0.01761040154983528
+-5.923662062192413e-5,0.01757158457833957
+-3.3613019432693456e-5,0.01755063790142304
+-2.1057361143545623e-5,0.017537845444160857
+-1.4102017705226928e-5,0.01752940755151119
+-9.919472449169899e-6,0.01752353136964117
+-7.247600005074752e-6,0.01751926793284237
+-5.458754315094613e-6,0.01751607316823318
+-4.2151006035299685e-6,0.01751361575179195
+-3.323232758355621e-6,0.01751168406986114
+-2.6667755990401974e-6,0.017510137616534138
+-2.1727568238066874e-6,0.01750888003257714
+-1.7937992750310765e-6,0.017507843397139875
+-1.4982125254477119e-6,0.017506978692230425
+-1.2642436477870115e-6,0.01750624980852978
+-1.0766187724705733e-6,0.01750562966201703
+-9.24396283448134e-7,0.01750509761099402
+-7.995968717103719e-7,0.01750463769771621
+-6.963076335628378e-7,0.017504237426426903
+-6.100832000933629e-7,0.017503886898331072
+-5.375374494886662e-7,0.01750357818896962
+-4.760601439147104e-7,0.017503304893260442
+-4.236170504643734e-7,0.017503061788459683
+-3.786068479257401e-7,0.017502844581325918
+-3.3975728204720806e-7,0.01750264971624999
+-3.0604885793705794e-7,0.004275364823468726

ChainOhmT/quadohmT.jl

@@ -1,13 +1,13 @@
 function quadfinT(N,i,uv)
 
     if (i>1 && i<mc)
-        xw = tr(uv,i)
+        xw = trr(uv,i)
         return xw
     end
 
     if mc==1
       if (AB[i,1]!=-Inf)&&(AB[i,2]!=Inf)
-        xw = tr(uv,i)
+        xw = trr(uv,i)
         return xw
       elseif AB[i,1]!=-Inf
           xw = rtr(uv,i)
@@ -21,7 +21,7 @@ function quadfinT(N,i,uv)
       end
     else
       if ((i==1 && AB[i,1]!=-Inf)||(i==mc && AB[i,2]!=Inf))
-        xw = tr(uv,i)
+        xw = trr(uv,i)
         return xw
       elseif i==1
         xw = ltr(uv,i)
@@ -33,7 +33,7 @@ function quadfinT(N,i,uv)
     return xw
 end
 
-function tr(t,i)
+function trr(t,i)
     s = Array{Float64}(undef, size(t))
     s[:,1] = ((AB[i,2]-AB[i,1])*t[:,1] .+ AB[i,2] .+ AB[i,1])./2
     s[:,2] = (AB[i,2]-AB[i,1]).*t[:,2].*wf(s[:,1],i)./2
@@ -65,9 +65,9 @@ function wf(x,i)
     if i==1
         y = 0
     elseif i==2
-        y = -pi*a*abs.(x) .* (coth.((beta/2).*x) .+ 1) #.* exp(-abs(x)/wc)
+        y = -( a*abs.(x).^s ./ wc^(s-1)) .* (coth.((beta/2).*x) .+ 1) #.* exp(-abs(x)/wc)
     elseif i==3
-        y = pi*a*abs.(x) .* (coth.((beta/2).*x) .+ 1) #.* exp(-abs(x)/wc)
+        y = ( a*abs.(x).^s ./ wc^(s-1)) .* (coth.((beta/2).*x) .+ 1) #.* exp(-abs(x)/wc)
     elseif i==4
         y = 0
     end

examples/puredephasing_temperature.jl

@@ -0,0 +1,110 @@
+#=
+    Example of a Pure Dephasing Model at zero 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 = 10 # 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(ω)
+
+β = 100
+
+curdir = @__DIR__
+
+dir_chaincoeff = abspath(joinpath(curdir, "../ChainOhmT/ohmicT"))
+
+#cpars is built here with the ChainOhmT routine for a thermalized ohmic 
+cpars  = readchaincoeffs("$dir_chaincoeff/chaincoeffs.h5",N,α,s,β) # 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 = 2 # 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 with 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 specified temperature 
+# (see The Theory of Open Quantum System, H.-P. Breuer & F. Petruccione 2002, Chapter 4)
+#----------
+
+Johmic(ω,s) = (2*α*ω^s)/(ωc^(s-1))
+
+time_analytical = LinRange(0.0,tfinal,Int(tfinal))
+
+Γohmic(t) = - quadgk(x -> Johmic(x,s)*(1 - cos(x*t))*coth(β*x/2)/x^2, 0, ωc)[1]
+
+Decoherence_ohmic(t) = 0.5*exp(Γohmic(t))
+
+#-------------
+# Plots
+#------------
+
+ρ12 = abs.(dat["data/Reduced ρ"][1,2,:])
+
+plot(time_analytical, t->Decoherence_ohmic(t), label="Analytics", title=L"Pure Dephasing, Ohmic $s=%$s$, $\beta = %$β ~\mathrm{K}$", linecolor=:black, xlabel="Time (arb. units)", ylabel=L"Coherence $|\rho_{12}(t)|$", linewidth=4, titlefontsize=16, legend=:best, legendfontsize=16, xguidefontsize=16, yguidefontsize=16, tickfontsize=10)
+
+plot!(dat["data/times"], ρ12, lw=4, ls=:dash, label="Numerics")
+

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, puredephasingmpo
+export chaincoeffs_ohmic, spinbosonmpo, methylbluempo, methylbluempo_correlated, methylbluempo_correlated_nocoupling, methylbluempo_nocoupling, ibmmpo, methylblue_S1_mpo, methylbluempo2, twobathspinmpo, xyzmpo, puredephasingmpo, tunnelingmpo
 
 export productstatemps, physdims, randmps, bonddims, elementmps