From daaf293d38381a7012880df3febf444f909ac6c5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Brieuc=20Le=20D=C3=A9?= Date: Thu, 25 Apr 2024 11:17:24 +0200 Subject: [PATCH] Write doc for puredephasing and time-dependency --- docs/make.jl | 2 +- docs/src/examples/puredephasing.md | 48 +++++++++++++++++------------- docs/src/examples/timedep.md | 11 ++++--- 3 files changed, 34 insertions(+), 27 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 7d483de..ac7b24a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -11,7 +11,7 @@ makedocs( pages = [ "index.md", "user-guide.md", - "Examples" => ["./examples/sbm.md", "./examples/puredephasing.md"], + "Examples" => ["./examples/sbm.md", "./examples/puredephasing.md", "./examples/timedep.md"], "theory.md", "Methods" => "methods.md", "dev.md" diff --git a/docs/src/examples/puredephasing.md b/docs/src/examples/puredephasing.md index 83480a4..42287f1 100644 --- a/docs/src/examples/puredephasing.md +++ b/docs/src/examples/puredephasing.md @@ -38,10 +38,12 @@ The T-TEDOPA method relies on a truncated chain mapping that transform the initi ## The code -First, we load the `MPSdynamics.jl` package to be able to perform the simulation, the `Plots.jl` one to plot the results, the `LaTeXStrings.jl` one to be able to use ``\LaTeX`` in the plots and eventually `QuadGK.jl` to perform the analytical integral calculations. +First, we load the `MPSdynamics.jl` package to be able to perform the simulation, the `Plots.jl` one to plot the results, the `LaTeXStrings.jl` one to be able to use ``\LaTeX`` in the plots and eventually `QuadGK.jl` to perform the analytical integral calculations. A notation is introduced for `+∞`. ```julia using MPSDynamics, Plots, LaTeXStrings, QuadGK + +const ∞ = Inf ``` We then define variables for the physical parameters of the simulation. Among these, two are convergence parameters: @@ -49,35 +51,43 @@ Among these, two 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 value of `β` determines whether the environment is thermalized or not. The example as it is is the zero-temperature case. For the finite-temperature case, 'β = ∞' has be commented instead of the line above that precises a `β` value. + ```julia #---------------------------- # Physical parameters #---------------------------- -ΔE = 0.008 # Energy of the electronic states - -d = 5 # number of Fock states of the chain modes +d = 10 # number of Fock states of the chain modes N = 30 # length of the chain α = 0.01 # coupling strength +ω0 = 0.008 # TLS gap + s = 1 # ohmicity ωc = 0.035 # Cut-off of the spectral density J(ω) -``` -The chain coefficient have still to be calculated. This part is the only difference in parameters for the zero-temperature case and the thermalized bath. For ``T = 0 ~ \text{K}``, the chain coeffficients can be calculated with the function [`MPSDynamics.chaincoeffs_ohmic`](@ref) : -```julia -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 -``` -Concerning the case ``T \neq 0 ~ \text{K}``, after having set up a temperature value, the chain coefficients can be calculated with the function [`MPSDynamics.chaincoeffs_finiteT`](@ref) +#β = 100 # Thermalized environment +β = ∞ # Case zero temperature T=0, β → ∞ +``` +The chain coefficient have then to be calculated. This part differs for the zero-temperature case and the thermalized bath. For ``β = ∞``, the chain coeffficients can be calculated with the function [`MPSDynamics.chaincoeffs_ohmic`](@ref) whereas the function [`MPSDynamics.chaincoeffs_finiteT`](@ref) is used when ``T \neq 0 ~ \text{K}``: ```julia -β = 100 - -cpars = chaincoeffs_finiteT(N, β; α=α, s=s, J=nothing, ωc=ωc, mc=4, mp=0, AB=nothing, iq=1, idelta=2, procedure=:Lanczos, Mmax=5000, save=false) # chain parameters, i.e. on-site energies ϵ_i, hopping energies t_i, and system-chain coupling c_0 +if β == ∞ + 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 +else + cpars = chaincoeffs_finiteT(N, β; α=α, s=s, J=nothing, ωc=ωc, mc=4, mp=0, AB=nothing, iq=1, idelta=2, procedure=:Lanczos, Mmax=5000, save=false) # chain parameters, i.e. on-site energies ϵ_i, hopping energies t_i, and system-chain coupling c_0 + #= #If cpars is stored in "../ChainOhmT/ohmicT" + curdir = @__DIR__ + dir_chaincoeff = abspath(joinpath(curdir, "../ChainOhmT/ohmicT")) + 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 + =# +end ``` +An option is provided if the coefficient are saved to be reused. 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)``. @@ -98,6 +108,8 @@ tfinal = 300.0 # simulation time method = :TDVP1 # time-evolution method +#method = :DTDVP # time-evolution method + D = 2 # MPS bond dimension ``` Using `MPSDynamics.jl` built-in methods we define the Pure Dephasing model MPO and the MPS representing the initial state. @@ -130,7 +142,7 @@ ob1 = OneSiteObservable("sz", sz, 1) # Simulation #------------ -A, dat = runsim(dt, tfinal, A, H; +A, dat = runsim(dt, tfinal, A, H, prec=1E-4; name = "pure dephasing model with temperature", method = method, obs = [ob1], @@ -143,7 +155,7 @@ A, dat = runsim(dt, tfinal, A, H; plot = true, ); ``` -After having performed the dynamics, the analytical decoherence function is numerically calculated with the help of the quadgk function (from the `QuadGK.jl` package). For the case ``T = 0 ~ \text{K}`` it reads +After having performed the dynamics, the analytical decoherence function is numerically calculated with the help of the quadgk function (from the `QuadGK.jl` package). It reads ```julia #---------- # Analytical results at specified temperature @@ -154,14 +166,10 @@ 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))/x^2, 0, ωc)[1] +Γ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)) ``` -whereas the analytical integral is changed for the case ``T \neq 0 ~ \text{K}`` -```julia -Γohmic(t) = - quadgk(x -> Johmic(x,s)*(1 - cos(x*t))*coth(β*x/2)/x^2, 0, ωc)[1] -``` Eventually, the stored reduced density matrix is compared against the analytical formula ```julia diff --git a/docs/src/examples/timedep.md b/docs/src/examples/timedep.md index 79b2885..17596ae 100644 --- a/docs/src/examples/timedep.md +++ b/docs/src/examples/timedep.md @@ -7,7 +7,7 @@ To simulate a drive or a laser pulse, a time-dependent Hamiltonian is often need ```math \hat{H}_\text{drive}(t) = \epsilon \hat{\sigma}_x \sin(\omega_\text{drive} t) ``` -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. +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. @@ -55,12 +55,11 @@ 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)``. -In this example we present two one-site implementation of TDVP that both preserves the unitarity of the evolution: +In this example we present only one-site implementation of TDVP that preserves the unitarity of the evolution: * the regular one-site method with the keyword `:TDVP1` where all the virtual bonds of the MPS have the same bond dimension ``D`` -* the adaptive method with the keyword `:DTDVP` where the bond dimension is locally increased at each time step if the TDVP projection error crosses a threshold value -Logically the constant bond dimension of the MPS for TDVP1 and the threshold of the projection error for DTDVP are their respective convergence parameter. +Logically the constant bond dimension of the MPS for TDVP1 is the respective convergence parameter. ```julia #----------------------- @@ -79,7 +78,7 @@ Using `MPSDynamics.jl` built-in methods we define the SBM MPO and the MPS repres 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. -In this part, the time-dependent terms of the MPO are stored in a list. This enables to add these terms to the MPO elements at each timestep and avoid the calculation of a new MPO. This way is not cumbersome since it adds a matrix sum at each time step. +In this part, the time-dependent terms of the MPO are stored in a list. This enables to add these terms to the MPO elements at each timestep and avoid the calculation of an entire new MPO. This way is not cumbersome since it adds a matrix sum at each time step. ```julia #--------------------------- @@ -97,7 +96,7 @@ H = spinbosonmpo(ω0, Δ, d, N, cpars) # MPO representation of the Hamiltonian A = productstatemps(physdims(H), state=[ψ, fill(unitcol(1,d), N)...]) # MPS representation of |ψ>|Vacuum> ``` -We then chose the observables that will be stored in the data and the [`MPSDynamics.runsim`](@ref) arguments. +We then choose the observables that will be stored in the data and the [`MPSDynamics.runsim`](@ref) arguments. ```julia #--------------------------- # Definition of observables