Skip to content

Commit a26d0c4

Browse files
BrieucLDBrieucLD
committed
Add Tailoring the Spectral Density example (+docs)
1 parent 517122c commit a26d0c4

File tree

4 files changed

+320
-2
lines changed

4 files changed

+320
-2
lines changed

‎docs/make.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ makedocs(
1212
"index.md",
1313
"user-guide.md",
1414
"nutshell.md",
15-
"Examples" => ["./examples/sbm.md", "./examples/puredephasing.md", "./examples/timedep.md", "./examples/anderson-model.md", "./examples/bath-observables.md", "./examples/protontransfer.md"],
15+
"Examples" => ["./examples/sbm.md", "./examples/puredephasing.md", "./examples/timedep.md", "./examples/anderson-model.md", "./examples/bath-observables.md", "./examples/protontransfer.md", "./examples/tailored_sd.md"],
1616
"convergence.md",
1717
"theory.md",
1818
"Methods" => "methods.md",

‎docs/src/examples/tailored_sd.md‎

Lines changed: 142 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
# Tailoring the Spectral Density
2+
3+
## Context
4+
5+
Here we detail an example to tailor the Spectral Density (SD) used to represent the bath. While the initial chain mapping procedure has been automatized for thermalized Ohmic type of SD (thanks to the methods [`MPSDynamics.chaincoeffs_ohmic`](@ref) and [`MPSDynamics.chaincoeffs_finiteT`](@ref)), it is possible to describe other types of SD. In this example, the SD is treated either with a continuous definition for several frequency domains or with a discrete description with frequency values and corresponding spectral density values as inputs.
6+
7+
The Spin-Boson Model is chosen to illustrate this example. For details about the model, we refer the interested reader to the example [The Spin-Boson Model](@ref).
8+
9+
## The code to define a specific spectral density
10+
11+
We describe here only differences with the Spin-Boson Model in order to parametrize the Hamiltonian with a specific spectral density.
12+
13+
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. A notation is introduced for `+∞`.
14+
15+
```julia
16+
using MPSDynamics, Plots, LaTeXStrings
17+
18+
const ∞ = Inf
19+
```
20+
We then define variables for the physical parameters of the simulation.
21+
22+
```julia
23+
#----------------------------
24+
# Physical parameters
25+
#----------------------------
26+
27+
ω0 = 0.2 # TLS gap
28+
29+
Δ = 0.0 # tunneling
30+
```
31+
32+
Now comes the definition of the spectral density. The type of SD has to be chosen between `:Discrete` and `:Continuous` (one choice has to be commented instead of the other).
33+
34+
For the `:Discrete` definition, the bath inputs become the discrete frequencies (`freqs`) with the corresponding spectral density values (`Jω`).
35+
36+
For the `:Continuous` definition, the bath inputs are as usual. Among these, two are convergence parameters:
37+
38+
* `d` is the number of states we retain for the truncated harmonic oscillators representation of environmental modes
39+
* `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
40+
41+
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.
42+
43+
After the parameters, the SD formulation has still to be defined according to the parameters. For that purpose, frequency domains have to be defined with the `AB` list and a function has to describe the spectral density formula for each domain. In order to follow the T-TEDOPA formulation, formula have to be modified as illustrated when the environment is thermalized. The presented example reproduces an Ohmic SD but frequency domains and formula can be changed at will.
44+
45+
```julia
46+
#-----------------------
47+
# Options and parameters for the Spectral Density
48+
# Definition of the frequency ranges and J(ω) formula for the continuous case
49+
#-----------------------
50+
51+
#Jω_type=:Discrete
52+
53+
Jω_type=:Continuous
54+
55+
if Jω_type==:Discrete
56+
57+
freqs = [0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50] # Frequency value
58+
59+
Jω = [0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10] # Value of the spectral density at the respective frequency
60+
61+
N = length(freqs) # The number of mode in the chain is set with the number of discrete frequencies
62+
63+
d = 6 # number of Fock states of the chain modes
64+
65+
elseif Jω_type==:Continuous
66+
67+
N = 30 # length of the chain
68+
69+
d = 6 # number of Fock states of the chain modes
70+
71+
α = 0.1 # coupling strength
72+
73+
s = 1 # ohmicity
74+
75+
ωc = 1.0 # Cut-off of the spectral density J(ω)
76+
77+
#β = 100 # Thermalized environment
78+
β = ∞ # Case zero temperature T=0, β → ∞
79+
80+
# AB segments correspond to the different frequency ranges where the spectral density is defined.
81+
# It corresponds to i=1 ; i=2 ; i=3 and i=4 of the Jω_fct(ω,i) function.
82+
83+
# The y formula and the frequency ranges can be changed while keeping the extra terms for thermalized environments.
84+
# This example reproduces an Ohmic type spectral density.
85+
86+
AB = [[-Inf -ωc];[-ωc 0];[0 ωc];[ωc Inf]]
87+
88+
function Jω_fct(ω,i)
89+
if i==1 # First segment of AB (here [-Inf -ωc])
90+
y = 0
91+
elseif i==2 # Second segment of AB (here [-ωc 0])
92+
if β == ∞
93+
y = 0 # zero for negative frequencies when β == Inf
94+
else
95+
y = -2*α*abs.(ω).^s / ωc^(s-1) .* (coth.((β/2).*ω) .+ 1) ./2 # .* (coth.((β/2).*ω) .+ 1) ./2 and the abs.(ω) have to be added when β != Inf (T-TEDOPA formulation)
96+
end
97+
elseif i==3 # Third segment of AB (here [0 ωc])
98+
if β == ∞
99+
y = 2*α*ω.^s / ωc^(s-1)
100+
else
101+
y = 2*α*ω.^s / ωc^(s-1) .* (coth.((β/2).*ω) .+ 1) ./2 # .* (coth.((β/2).*ω) .+ 1) ./2 has to be added when β != Inf (T-TEDOPA formulation)
102+
end
103+
elseif i==4 # Fourth segment of AB (here [ωc Inf])
104+
y = 0
105+
end
106+
return y
107+
end
108+
109+
else
110+
111+
throw(ErrorException("The spectral density has either to be Continuous or Discrete"))
112+
113+
end
114+
```
115+
116+
Now that the type of SD and the corresponding parameters have been filled in, the chain parameters have to be calculated.
117+
118+
For both cases, the method [`MPSDynamics.chaincoeffs_finiteT`](@ref) can be called with the respective inputs.
119+
120+
```julia
121+
#-----------------------
122+
# Calculation of chain parameters
123+
#-----------------------
124+
125+
if Jω_type==:Discrete
126+
127+
cpars = chaincoeffs_finiteT_discrete(β, freqs, Jω; procedure=:Lanczos, Mmax=5000, save=false) # chain parameters, i.e. on-site energies ϵ_i, hopping energies t_i, and system-chain coupling c_0
128+
129+
#= If cpars is stored in "../ChainOhmT/discreteT"
130+
curdir = @__DIR__
131+
dir_chaincoeff = abspath(joinpath(curdir, "../ChainOhmT/discreteT"))
132+
cpars = readchaincoeffs("$dir_chaincoeff/chaincoeffs.h5", N, β) # chain parameters, i.e. on-site energies ϵ_i, hopping energies t_i, and system-chain coupling c_0
133+
=#
134+
135+
elseif Jω_type==:Continuous
136+
137+
cpars = chaincoeffs_finiteT(N, β, false; α=α, s=s, J=Jω_fct, ωc=ωc, mc=size(AB)[1], mp=0, AB=AB, 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
138+
139+
end
140+
```
141+
142+
Next, similarly to the spin-boson example, the simulation parameters are set with the initial MPO and MPS in order to carry out the dynamics propagation. Therefore, this modified example shows an easy way to tailor the spectral density for any system and Hamiltonian.

‎docs/src/index.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ The package may be installed by typing `]` in a Julia REPL to enter the package
1616
## Table of Contents
1717

1818
```@contents
19-
Pages = ["index.md", "user-guide.md", "examples/sbm.md", "examples/puredephasing.md", "examples/timedep.md", "examples/anderson-model.md", "examples/bath-observables.md", "examples/protontransfer.md", "convergence.md", "theory.md", "methods.md", "dev.md"]
19+
Pages = ["index.md", "user-guide.md", "examples/sbm.md", "examples/puredephasing.md", "examples/timedep.md", "examples/anderson-model.md", "examples/bath-observables.md", "examples/protontransfer.md", "examples/tailored_sd.md", "convergence.md", "theory.md", "methods.md", "dev.md"]
2020
Depth = 3
2121
```
2222

‎examples/sbm_tailored_specdensity.jl‎

Lines changed: 176 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,176 @@
1+
#=
2+
Example of a zero-temperature Spin-Boson Model with a tailored spectral density J(ω). This spectral density can be either defined with discrete
3+
modes or continuous formula.
4+
5+
The dynamics is simulated using the T-TEDOPA method that maps the normal modes environment into a non-uniform tight-binding chain.
6+
7+
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
8+
9+
Two variants of the one-site Time Dependent Variational Principal (TDVP) are presented for the time evolution of the quantum state.
10+
=#
11+
12+
using MPSDynamics, Plots, LaTeXStrings
13+
14+
const ∞ = Inf
15+
#----------------------------
16+
# Physical parameters
17+
#----------------------------
18+
19+
ω0 = 0.2 # TLS gap
20+
21+
Δ = 0.0 # tunneling
22+
23+
#-----------------------
24+
# Options and parameters for the Spectral Density
25+
# Definition of the frequency ranges and J(ω) formula for the continuous case
26+
#-----------------------
27+
28+
#Jω_type=:Discrete
29+
30+
Jω_type=:Continuous
31+
32+
if Jω_type==:Discrete
33+
34+
freqs = [0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50] # Frequency value
35+
36+
Jω = [0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10] # Value of the spectral density at the respective frequency
37+
38+
N = length(freqs) # The number of mode in the chain is set with the number of discrete frequencies
39+
40+
d = 6 # number of Fock states of the chain modes
41+
42+
elseif Jω_type==:Continuous
43+
44+
N = 30 # length of the chain
45+
46+
d = 6 # number of Fock states of the chain modes
47+
48+
α = 0.1 # coupling strength
49+
50+
s = 1 # ohmicity
51+
52+
ωc = 1.0 # Cut-off of the spectral density J(ω)
53+
54+
#β = 100 # Thermalized environment
55+
β = ∞ # Case zero temperature T=0, β → ∞
56+
57+
# AB segments correspond to the different frequency ranges where the spectral density is defined.
58+
# It corresponds to i=1 ; i=2 ; i=3 and i=4 of the Jω_fct(ω,i) function.
59+
60+
# The y formula and the frequency ranges can be changed while keeping the extra terms for thermalized environments.
61+
# This example reproduces an Ohmic type spectral density.
62+
63+
AB = [[-Inf -ωc];[-ωc 0];[0 ωc];[ωc Inf]]
64+
65+
function Jω_fct(ω,i)
66+
if i==1 # First segment of AB (here [-Inf -ωc])
67+
y = 0
68+
elseif i==2 # Second segment of AB (here [-ωc 0])
69+
if β == ∞
70+
y = 0 # zero for negative frequencies when β == Inf
71+
else
72+
y = -2*α*abs.(ω).^s / ωc^(s-1) .* (coth.((β/2).*ω) .+ 1) ./2 # .* (coth.((β/2).*ω) .+ 1) ./2 and the abs.(ω) have to be added when β != Inf (T-TEDOPA formulation)
73+
end
74+
elseif i==3 # Third segment of AB (here [0 ωc])
75+
if β == ∞
76+
y = 2*α*ω.^s / ωc^(s-1)
77+
else
78+
y = 2*α*ω.^s / ωc^(s-1) .* (coth.((β/2).*ω) .+ 1) ./2 # .* (coth.((β/2).*ω) .+ 1) ./2 has to be added when β != Inf (T-TEDOPA formulation)
79+
end
80+
elseif i==4 # Fourth segment of AB (here [ωc Inf])
81+
y = 0
82+
end
83+
return y
84+
end
85+
86+
else
87+
88+
throw(ErrorException("The spectral density has either to be Continuous or Discrete"))
89+
90+
end
91+
92+
#-----------------------
93+
# Calculation of chain parameters
94+
#-----------------------
95+
96+
if Jω_type==:Discrete
97+
98+
cpars = chaincoeffs_finiteT_discrete(β, freqs, Jω; procedure=:Lanczos, Mmax=5000, save=false) # chain parameters, i.e. on-site energies ϵ_i, hopping energies t_i, and system-chain coupling c_0
99+
100+
#= If cpars is stored in "../ChainOhmT/discreteT"
101+
curdir = @__DIR__
102+
dir_chaincoeff = abspath(joinpath(curdir, "../ChainOhmT/discreteT"))
103+
cpars = readchaincoeffs("$dir_chaincoeff/chaincoeffs.h5", N, β) # chain parameters, i.e. on-site energies ϵ_i, hopping energies t_i, and system-chain coupling c_0
104+
=#
105+
106+
elseif Jω_type==:Continuous
107+
108+
cpars = chaincoeffs_finiteT(N, β, false; α=α, s=s, J=Jω_fct, ωc=ωc, mc=size(AB)[1], mp=0, AB=AB, 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
109+
110+
end
111+
112+
113+
#-----------------------
114+
# Simulation parameters
115+
#-----------------------
116+
117+
dt = 0.5 # time step
118+
119+
tfinal = 30.0 # simulation time
120+
121+
method = :TDVP1 # Regular one-site TDVP (fixed bond dimension)
122+
123+
# method = :DTDVP # Adaptive one-site TDVP (dynamically updating bond dimension)
124+
125+
convparams = [2,4,6] # MPS bond dimension (1TDVP)
126+
127+
# convparams = [1e-2, 1e-3, 1e-4] # threshold value of the projection error (DTDVP)
128+
129+
#---------------------------
130+
# MPO and initial state MPS
131+
#---------------------------
132+
133+
H = spinbosonmpo(ω0, Δ, d, N, cpars) # MPO representation of the Hamiltonian
134+
135+
ψ = unitcol(1,2) # Initial up-z system state
136+
137+
A = productstatemps(physdims(H), state=[ψ, fill(unitcol(1,d), N)...]) # MPS representation of |ψ>|Vacuum>
138+
139+
#---------------------------
140+
# Definition of observables
141+
#---------------------------
142+
143+
ob1 = OneSiteObservable("sz", sz, 1)
144+
145+
ob2 = OneSiteObservable("chain mode occupation", numb(d), (2,N+1))
146+
147+
ob3 = TwoSiteObservable("SXdisp", sx, MPSDynamics.disp(d), [1], collect(2:N+1))
148+
149+
#-------------
150+
# Simulation
151+
#------------
152+
153+
A, dat = runsim(dt, tfinal, A, H;
154+
name = "ohmic spin boson model",
155+
method = method,
156+
obs = [ob2,ob3],
157+
convobs = [ob1],
158+
params = @LogParams(Δ, ω0, N, Jω_type),
159+
convparams = convparams,
160+
verbose = false,
161+
savebonddims = true, # this keyword argument enables the bond dimension at each time step to be saved when using DTDVP
162+
save = true,
163+
plot = true,
164+
);
165+
166+
#----------
167+
# Plots
168+
#----------
169+
170+
display(method == :TDVP1 && plot(dat["data/times"], dat["convdata/sz"], label=["Dmax = 2" "Dmax = 4" "Dmax = 6"], xlabel=L"t",ylabel=L"\sigma_z"))
171+
172+
method == :DTDVP && plot(dat["data/times"], dat["convdata/sz"], label=["p = 1e-2" "p = 1e-3" "p = 1e-4"], xlabel=L"t",ylabel=L"\sigma_z")
173+
174+
method == :DTDVP && heatmap(dat["data/times"], collect(0:N+1), dat["data/bonddims"], xlabel=L"t",ylabel="bond index")
175+
176+
heatmap(dat["data/times"], collect(1:N), abs.(dat["data/SXdisp"][1,:,:]), xlabel=L"t",ylabel="chain mode")

0 commit comments

Comments
 (0)