From 78f952b07039bfc34f43e6afc0237f8367e64697 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Brieuc=20Le=20D=C3=A9?= Date: Wed, 12 Jun 2024 16:05:45 +0200 Subject: [PATCH] Fix typos and add docs --- docs/src/examples/anderson-model.md | 19 +-- docs/src/examples/bath-observables.md | 6 +- src/deprecated/fundamentals_deprecated.jl | 16 +++ src/deprecated/models_deprecated.jl | 18 +++ src/fundamentals.jl | 165 +++++++++++++--------- src/models.jl | 75 ++++------ src/mpsBasics.jl | 4 +- 7 files changed, 173 insertions(+), 130 deletions(-) create mode 100644 src/deprecated/fundamentals_deprecated.jl create mode 100644 src/deprecated/models_deprecated.jl diff --git a/docs/src/examples/anderson-model.md b/docs/src/examples/anderson-model.md index db42a96..7c99ac8 100644 --- a/docs/src/examples/anderson-model.md +++ b/docs/src/examples/anderson-model.md @@ -1,6 +1,6 @@ # The Anderson Impurity Model -Here we give some context on the examples provided in `MPSDynamics/examples/anderson_model_double.jl` and `MPSDynamics/examples/anderson_model_interleaved.jl`. In these two examples, we use the fermionic chain mapping proposed in [^khon_efficient_2021] to show how to perform tensor network simulations of the Single Impurity Anderson Model (SIAM) with the `MPSDynamics.jl` library. +Here we give some context on the examples provided in `MPSDynamics/examples/anderson_model_double.jl` and `MPSDynamics/examples/anderson_model_interleaved.jl`. In these two examples, we use the fermionic chain mapping proposed in [^khon_efficient_2021] to show how to perform tensor network simulations of the [Single Impurity Anderson Model](https://en.wikipedia.org/wiki/Anderson_impurity_model) (SIAM) with the `MPSDynamics.jl` library. ## Basics of the fermionic chain mapping @@ -8,9 +8,10 @@ Here we give some context on the examples provided in `MPSDynamics/examples/ande Before giving the code implementation, a short recap on the problem statement and on the fermionic mapping used to solved it. The SIAM Hamiltonian is defined as: ```math - \hat H^\text{SIAM} = \hat H_\text{loc} + \hat H_\text{hyb} + \hat H_\text{cond} = \overbrace{\epsilon_d \hat d^\dagger \hat d}^{\hat H_\text{loc}} + \underbrace{\sum_{k} V_k \Big( \hat d^\dagger \hat c_k + \hat c_k^\dagger \hat d \Big)}_{H_\text{hyb}} + \underbrace{\sum_k \epsilon_k \hat c_k^\dagger \hat c_k}_{H_I^\text{chain}}. + \hat H^\text{SIAM} = \hat H_\text{loc} + \hat H_\text{hyb} + \hat H_\text{cond} = \overbrace{\epsilon_d \hat d^\dagger \hat d}^{\hat H_\text{loc}} + \underbrace{\sum_{k} V_k \Big( \hat d^\dagger \hat c_k + \hat c_k^\dagger \hat d \Big)}_{\hat H_\text{hyb}} + \underbrace{\sum_k \epsilon_k \hat c_k^\dagger \hat c_k}_{\hat H_\text{cond}}. ``` -All of the operators obey to the usual fermionic anti-commutation relations: $\{\hat c_i, \hat c_j^\dagger \} = \delta_{ij}$, $\{\hat c_i, \hat c_j \} =\{\hat c_i^\dagger, \hat c_j^\dagger \} =0$ $\forall i,j$. The chain mapping is based on a thermofield-like transformation + +This Hamiltonian represents a located impurity ($\hat H_\text{loc}$), conduction electrons ($\hat H_\text{cond}$) and a hybridization term between the impurity and the conduction electrons ($\hat H_\text{hyb}$). All of the operators obey to the usual fermionic anti-commutation relations: $\{\hat c_i, \hat c_j^\dagger \} = \delta_{ij}$, $\{\hat c_i, \hat c_j \} =\{\hat c_i^\dagger, \hat c_j^\dagger \} =0$ $\forall i,j$. The chain mapping is based on a thermofield-like transformation [^devega_thermo_2015], performed with fermions: ancillary fermionic operators $\hat c_{2k}$ are defined, one for each of the original fermionic modes $\hat c_{1k}$. A Bogoliubov transformation is then applied, so that two new fermionic modes $\hat f_{1k}$ and $\hat f_{2k}$ are defined as a linear combination of $\hat c_{1k}$ and $\hat c_{2k}$. Two chains are defined: the chain labelled $1$ for the empty modes: ```math \hat f_{1k}=e^{-iG} \hat c_k e^{iG}= \cosh(\theta_k) \hat c_{1k} -\sinh(\theta_k) \hat c_{2k}^\dagger \\ @@ -136,7 +137,7 @@ where the matrices are defined as: ## Code implementation -The fermionic mapping can be strightforwardly implemented with the methods of `MPSDyanmics.jl`. We provide two examples: +The fermionic mapping can be strightforwardly implemented with the methods of `MPSDynamics.jl`. We provide two examples: - `examples/anderson_model_double`, simulating the SIAM with the double-chain geometry - `examples/anderson_model_interleaved`, simulating the SIAM with the interleaved-chain geometry They only differ on the way the Hamiltonian MPO is defined. We now briefly review the important bits of the code. @@ -149,7 +150,7 @@ N = 40 # number of chain sites Ed = 0.3 # energy of the impurity ϵd = Ed - μ # energy of the impurity minus the chemical potential ``` -The `chaincoeffs_fermionic` function is needed to compute the chain coefficients. It requires as inputs the number of modes of each chain `N`, the inverse temperature `β`, a label to specify if the chain modes are empty (label is `1.0`) or filled (label is `2.0`), and both the dispersion relation $\epsilon_k$ and the fermionic spectral density funciton $V_k$. +The [`MPSDynamics.chaincoeffs_fermionic`](@ref) function is needed to compute the chain coefficients. It requires as inputs the number of modes of each chain `N`, the inverse temperature `β`, a label to specify if the chain modes are empty (label is `1.0`) or filled (label is `2.0`), and both the dispersion relation $\epsilon_k$ and the fermionic spectral density funciton $V_k$. ```julia function ϵ(x) return x @@ -174,7 +175,7 @@ prec = 0.0001 # precision for the adaptive TDVP and with this we are ready to construct the Hamiltonian MPO and specify the initial state, which will obviously differ depending on the chosen geometry. #### Double chain geometry -The Hamiltonian is defined using the `tightbinding_mpo` function, which takes as an input the number of modes of each chain `N`, the defect's energy `ϵd`, and the chain coefficients of the first `chainparams1` and second `chainparams2` chain. The MPS for the initial state is a factorized state made of: N filled states, a filled impurity, and N empty states. +The Hamiltonian is defined using the [`MPSDynamics.tightbinding_mpo`](@ref) function, which takes as an input the number of modes of each chain `N`, the defect's energy `ϵd`, and the chain coefficients of the first `chainparams1` and second `chainparams2` chain. The MPS for the initial state is a factorized state made of: N filled states, a filled impurity, and N empty states. ```julia H = tightbinding_mpo(N, ϵd, chainparams1, chainparams2) @@ -198,7 +199,7 @@ A, dat = runsim(dt, T, A, H; method = method, obs = [ob1, ob2, ob3], convobs = [ob1], - params = @LogParams(N, ϵd, β, c1, c2), + params = @LogParams(N, ϵd, β), convparams = [prec], Dlim = Dmax, savebonddims = true, # we want to save the bond dimension @@ -276,7 +277,7 @@ plot(p2, p3, p4, p5, p1, layout = (3, 2), size = (1400, 1200)) #### Interleaved chain geometry -The Hamiltonian is defined using the `interleaved_tightbinding_mpo` function, which takes as an input the number of modes of each chain `N`, the defect's energy `ϵd`, and the chain coefficients of the first `chainparams1` and second `chainparams2` chain. The MPS for the initial state is a factorized state (bond dimension 1) made of: a filled impurity, and 2N alternate filled-empty states. +The Hamiltonian is defined using the [`MPSDynamics.interleaved_tightbinding_mpo`](@ref) function, which takes as an input the number of modes of each chain `N`, the defect's energy `ϵd`, and the chain coefficients of the first `chainparams1` and second `chainparams2` chain. The MPS for the initial state is a factorized state (bond dimension 1) made of: a filled impurity, and 2N alternate filled-empty states. ```julia H = interleaved_tightbinding_mpo(N, ϵd, chainparams1, chainparams2) @@ -306,7 +307,7 @@ A, dat = runsim(dt, T, A, H; method = method, obs = [ob1, ob2], convobs = [ob1], - params = @LogParams(N, ϵd, β, c1, c2), + params = @LogParams(N, ϵd, β), convparams = [prec], Dlim = Dmax, savebonddims = true, # we want to save the bond dimension diff --git a/docs/src/examples/bath-observables.md b/docs/src/examples/bath-observables.md index 050682c..4ba0904 100644 --- a/docs/src/examples/bath-observables.md +++ b/docs/src/examples/bath-observables.md @@ -79,7 +79,7 @@ for t in 1:T end ``` -During a numerical simulation we work with chain of finite length $N$. This truncation on chain modes, introduces a sampling on the modes in the original star-like environment. To recover the frequency modes that are implicitly sampled, one has to diagonalize the tri-diagonal $N\times N$ matrix $H^\text{chain}$, where the diagonal is formed by the $e_n$ coefficients, that is the chain's frequencies, and the upper and lower diagonals by the $N-1$ hopping coefficients $t_n$. The unitary matrix that defines the change of basis from the star-like to the chain-like environment is $U_n$, constituted by the eigenvectors of $H^\text{chain}$. In the code, we use the `eigenchain` function: +During a numerical simulation we work with chain of finite length $N$. This truncation on chain modes, introduces a sampling on the modes in the original star-like environment. To recover the frequency modes that are implicitly sampled, one has to diagonalize the tri-diagonal $N\times N$ matrix $H^\text{chain}$, where the diagonal is formed by the $e_n$ coefficients, that is the chain's frequencies, and the upper and lower diagonals by the $N-1$ hopping coefficients $t_n$. The unitary matrix that defines the change of basis from the star-like to the chain-like environment is $U_n$, constituted by the eigenvectors of $H^\text{chain}$. In the code, we use the [`MPSDynamics.eigenchain`](@ref) function: ```julia omeg = eigenchain(cpars, nummodes=N).values @@ -94,7 +94,7 @@ Therefore to calculate the number of modes of the environment associated to a sp \end{aligned} ``` -This is done in the code using the `measuremodes(X, cpars[1], cpars[2])` function, which outputs the vector of the diagonal elements of the operators, in the following way: +This is done in the code using the [`MPSDynamics.measuremodes`](@ref) function, which outputs the vector of the diagonal elements of the operators, in the following way: ```julia bath_occup = mapslices(X -> measuremodes(X, cpars[1], cpars[2]), dat["data/cdagc"], dims = [1,2]) @@ -126,7 +126,7 @@ It is possible to invert the thermofield transformation (details in [^riva_therm \end{aligned} ``` -We remark that in the thermofield case, a negative frequency $\omega_{2k}$ is associated to each positive frequency $\omega_{1k}$. The sampling is therefore symmetric around zero. This marks a difference with T-TEDOPA, where the sampling of frequencies was obtained through the thermalized measure $d\mu(\beta) = \sqrt{J(\omega, \beta)}d\omega$, and was not symmetric. To recover the results for the physical bath of frequencies starting from the results of our simulations, that were conducted using the T-TEDOPA chain mapping, we need to do an extrapolation for all of the mean values, in order to have their values for each $\omega$ at $-\omega$ as well. This is done in the code with the `physical_occup` function: +We remark that in the thermofield case, a negative frequency $\omega_{2k}$ is associated to each positive frequency $\omega_{1k}$. The sampling is therefore symmetric around zero. This marks a difference with T-TEDOPA, where the sampling of frequencies was obtained through the thermalized measure $d\mu(\beta) = \sqrt{J(\omega, \beta)}d\omega$, and was not symmetric. To recover the results for the physical bath of frequencies starting from the results of our simulations, that were conducted using the T-TEDOPA chain mapping, we need to do an extrapolation for all of the mean values, in order to have their values for each $\omega$ at $-\omega$ as well. This is done in the code with the [`MPSDynamics.physical_occup`](@ref) function: ```julia bath_occup_phys = physical_occup(correlations_cdag[:,:,T], correlations_c[:,:,T], omeg, bath_occup[:,:,T], β, N) diff --git a/src/deprecated/fundamentals_deprecated.jl b/src/deprecated/fundamentals_deprecated.jl new file mode 100644 index 0000000..46faa59 --- /dev/null +++ b/src/deprecated/fundamentals_deprecated.jl @@ -0,0 +1,16 @@ +# Doesn't really work because indices change order in ITensors +function MPStoVector(mps::MPS) + N = length(mps) + A = [Array(mps[i], mps[i].inds...) for i=1:N] + dims = size(A[1]) + A[1] = reshape(A[1], 1, dims...) + A[1] = permutedims(A[1], [1,3,2]) + for i=2:N-1 + A[i] = permutedims(A[i], [3,1,2]) + end + dims = size(A[N]) + A[N] = reshape(A[N], dims..., 1) + A[N] = permutedims(A[N], [1,3,2]) + return A +end + diff --git a/src/deprecated/models_deprecated.jl b/src/deprecated/models_deprecated.jl new file mode 100644 index 0000000..dada732 --- /dev/null +++ b/src/deprecated/models_deprecated.jl @@ -0,0 +1,18 @@ +#Deprecated# +function getchaincoeffs(nummodes, α, s, beta, ωc=1) + matlabdir = ENV["MATDIR"] + astr = paramstring(α, 2) + bstr = paramstring(beta, 3) + datfname = string("jacerg_","a",astr,"s$s","beta",bstr,".csv") + chaincoeffs = readdlm(string(matlabdir,datfname),',',Float64,'\n') + es = ωc .* chaincoeffs[:,1] + ts = ωc .* chaincoeffs[1:end-1,2] + c0 = ωc * sqrt(chaincoeffs[end,2]/pi) + Nmax = length(es) + if Nmax < nummodes + throw(ErrorException("no data for nummodes > $Nmax")) + end + return [es[1:nummodes], ts[1:nummodes-1], c0] +end +## + diff --git a/src/fundamentals.jl b/src/fundamentals.jl index 53ddef5..bf82bc6 100644 --- a/src/fundamentals.jl +++ b/src/fundamentals.jl @@ -4,7 +4,7 @@ crea(d) = diagm(-1 => [sqrt(i) for i=1:d-1]) anih(d) = Matrix(crea(d)') numb(d) = crea(d)*anih(d) """ - function disp(d) + disp(d) Mass and frequency-weighted displacement operator `` @@ -14,13 +14,12 @@ X = \\frac{1}{2}(a + a^{\\dagger}) """ disp(d) = (1/sqrt(2))*(crea(d)+anih(d)) """ - function disp(d,ωvib,m) + disp(d,ωvib,m) Displacement operator `` X = \\frac{\\sqrt{2}}{2\\sqrt{m \\omega_{vib}}}(a + a^{\\dagger}) `` - """ disp(d,ωvib,m) = (1/(2*sqrt(m*ωvib/2)))*(crea(d)+anih(d)) mome(d) = (1/sqrt(2))*(im*(crea(d)-anih(d))) @@ -93,7 +92,7 @@ function choose(x, n) end """ - function therHam(psi, site1, site2) + therHam(psi, site1, site2) Calculates Hβ such that ρ = e^(-βH) for some density matrix ρ obatined from tracing out everything outside the range [site1,site2] in the MPS psi """ @@ -139,6 +138,10 @@ function savecsv(dir, fname, dat) writedlm(string(dir, fname,".csv"), dat, ',') end +""" + eigenchain(cparams; nummodes=nothing) + +""" function eigenchain(cparams; nummodes=nothing) if nummodes==nothing nummodes = length(cparams[1]) @@ -149,6 +152,10 @@ function eigenchain(cparams; nummodes=nothing) return eigen(hmat) end +""" + thermaloccupations(β, cparams...) + +""" function thermaloccupations(β, cparams...) es=cparams[1] ts=cparams[2] @@ -159,11 +166,15 @@ function thermaloccupations(β, cparams...) ((U)^2) * (1 ./ (exp.(β.*S).-1)) end +""" + measuremodes(A, chainsection::Tuple{Int64,Int64}, e::Array{Float64,1}, t::Array{Float64,1}) + +""" function measuremodes(A, chainsection::Tuple{Int64,Int64}, e::Array{Float64,1}, t::Array{Float64,1}) N=abs(chainsection[1]-chainsection[2])+1 e=e[1:N] t=t[1:N-1] - d=size(A[chainsection[1]])[2]#assumes constant number of Foch states + d=size(A[chainsection[1]])[2]#assumes constant number of Fock states bd=crea(d) b=anih(d) hmat = diagm(0=>e, 1=>t, -1=>t) @@ -171,13 +182,26 @@ function measuremodes(A, chainsection::Tuple{Int64,Int64}, e::Array{Float64,1}, U = F.vectors return real.(diag(U' * measure2siteoperator(A, bd, b, chainsection, conj=true, herm_cis=true) * U)) end -#for longer chains it can be worth calculating U in advance + +""" + measuremodes(A, chainsection::Tuple{Int64,Int64}, U::AbstractArray) + +for longer chains it can be worth calculating U in advance + + +""" function measuremodes(A, chainsection::Tuple{Int64,Int64}, U::AbstractArray) - d=size(A[chainsection[1]])[2]#assumes constant number of Foch states + d=size(A[chainsection[1]])[2]#assumes constant number of Fock states bd=crea(d) b=anih(d) return real.(diag(U' * measure2siteoperator(A, bd, b, chainsection, conj=true, herm_cis=true) * U)) end + +""" + measuremodes(adaga, e=1.0, t=1.0) + + +""" function measuremodes(adaga, e=1.0, t=1.0) N = size(adaga)[1] hmat = diagm(0=>fill(e,N), 1=>fill(t,N-1), -1=>fill(t,N-1)) @@ -185,6 +209,12 @@ function measuremodes(adaga, e=1.0, t=1.0) U = F.vectors return real.(diag(U' * adaga * U)) end + +""" + measuremodes(adaga, e::Vector, t::Vector) + + +""" function measuremodes(adaga, e::Vector, t::Vector) N = size(adaga)[1] hmat = diagm(0=>e[1:N], 1=>t[1:N-1], -1=>t[1:N-1]) @@ -196,24 +226,24 @@ end """ measurecorrs(oper, , e::Vector, t::Vector) - ### Parameters +### Parameters - `oper``: Square matrix (Matrix{Float64}) representing the operator to be transformed. - `e``: Vector (Vector{Float64}) of diagonal (on-site energy) chain coefficients. - `t`: Vector (Vector{Float64}) of off-diagonal (hopping terms) chain coefficients. +`oper`: Square matrix (Matrix{Float64}) representing the operator to be transformed. +`e`: Vector (Vector{Float64}) of diagonal (on-site energy) chain coefficients. +`t`: Vector (Vector{Float64}) of off-diagonal (hopping terms) chain coefficients. - ### Returns +### Returns - Matrix{Float64}: This matrix is the operator `oper` transformed back from the chain - representation to the representation corresponding to the extended bath. The resulting - operator represents quantities like mode occupations or other properties in the basis - of environmental modes associated with specific frequencies ``\\omega_i``. +Matrix{Float64}: This matrix is the operator `oper` transformed back from the chain +representation to the representation corresponding to the extended bath. The resulting +operator represents quantities like mode occupations or other properties in the basis +of environmental modes associated with specific frequencies ``\\omega_i``. - ### Description +### Description - This function performs a basis transformation of the operator `oper`. Specifically, - this transformation reverses the unitary transformation that maps the extended bath - Hamiltonian into the chain representation. +This function performs a basis transformation of the operator `oper`. Specifically, +this transformation reverses the unitary transformation that maps the extended bath +Hamiltonian into the chain representation. """ function measurecorrs(oper, e::Vector, t::Vector) @@ -228,15 +258,15 @@ end """ cosineh(omega, bet) - Calculates the hyperbolic cosine function function based on the input parameters, - for the Bogoliubov transformation necessary for the thermofield transformation. +Calculates the hyperbolic cosine function function based on the input parameters, +for the Bogoliubov transformation necessary for the thermofield transformation. - # Arguments - - `omega::Float64`: The frequency parameter. - - `bet::Float64`: The beta parameter. +# Arguments +- `omega::Float64`: The frequency parameter. +- `bet::Float64`: The beta parameter. - # Returns - - `Float64`: The result of the modified cosine function. +# Returns +- `Float64`: The result of the modified cosine function. """ function cosineh(omega, bet) return 1/sqrt(1 - exp(-omega * (bet))) @@ -245,15 +275,15 @@ end """ sineh(omega, bet) - Calculates the hyperbolic sine function function based on the input parameters, - for the Bogoliubov transformation necessary for the thermofield transformation. +Calculates the hyperbolic sine function function based on the input parameters, +for the Bogoliubov transformation necessary for the thermofield transformation. - # Arguments - - `omega::Float64`: The frequency parameter. - - `bet::Float64`: The beta parameter. +# Arguments +- `omega::Float64`: The frequency parameter. +- `bet::Float64`: The beta parameter. - # Returns - - `Float64`: The result of the modified cosine function. +# Returns +- `Float64`: The result of the modified cosine function. """ function sineh(omega, bet) return 1/sqrt(-1 + exp(omega * float(bet))) @@ -262,20 +292,20 @@ end """ physical_occup(corr_constr, corr_destr, omega, occup, b, M) - Calculates the physical occupation based on correlation matrices, omega values, - and other parameters. The physical occupation in the original frequency environment - is computed by reverting the thermofield transformation. +Calculates the physical occupation based on correlation matrices, omega values, +and other parameters. The physical occupation in the original frequency environment +is computed by reverting the thermofield transformation. - # Arguments - - `corr_constr::Matrix{ComplexF64}`: The correlation construction matrix. - - `corr_destr::Matrix{ComplexF64}`: The correlation destruction matrix. - - `omega::Vector{Float64}`: The omega values. - - `occup::Matrix{Float64}`: The occupation matrix. - - `b::Float64`: The beta parameter. - - `M::Int`: The number of points for interpolation. +# Arguments +- `corr_constr::Matrix{ComplexF64}`: The correlation construction matrix. +- `corr_destr::Matrix{ComplexF64}`: The correlation destruction matrix. +- `omega::Vector{Float64}`: The omega values. +- `occup::Matrix{Float64}`: The occupation matrix. +- `b::Float64`: The beta parameter. +- `M::Int`: The number of points for interpolation. - # Returns - - `Vector{Float64}`: The physical occupation values. +# Returns +- `Vector{Float64}`: The physical occupation values. """ function physical_occup(corr_constr, corr_destr, omega, occup, b, M) x = range(-1, stop=1, length=M) @@ -309,7 +339,7 @@ end """ - findchainlength(T, cparams...; eps=10^-6) + findchainlength(T, cparams; eps=10^-6, verbose=false) Estimate length of chain required for a particular set of chain parameters by calculating how long an excitation on the first site takes to reach the end. The chain length is given as the length required for the excitation to have just @@ -339,7 +369,7 @@ function findchainlength(T, cparams; eps=10^-4, verbose=false) end """ - chainprop(t, cparams...) + chainprop(t, cparams) Propagate an excitation placed initially on the first site of a tight-binding chain with parameters given by cparams for a time t and return occupation expectation for each site. @@ -355,6 +385,11 @@ function chainprop(t, cparams) [real.(transpose(U[:,1].*exp.(im*t.*S))*U[:,i]*transpose(U[:,i])*(U[:,1].*exp.(-im*t.*S))) for i in 1:N] end +""" + endsiteocc(t, cparams) + + +""" function endsiteocc(t, cparams) es=cparams[1] ts=cparams[2] @@ -389,7 +424,7 @@ end sd(nums...) = sd([nums...]) """ - rmsd(dat1::Vector{Float64}, dat2::Vector{Float64}) + rmsd(ob1, ob2) Calculate the root mean squared difference between two measurements of an observable over the same time period. @@ -407,12 +442,13 @@ end Calculate complete dynamical map to time step at which ps1, ps2, ps3 and ps4 are specified. +# Arguments +- `ps1` : time evolved system density matrix starting from initial state up +- `ps2` : time evolved system density matrix starting from initial state down +- `ps3` : time evolved system density matrix starting from initial state (up + down)/sqrt(2) +- `ps4` : time evolved system density matrix starting from initial state (up - i*down)/sqrt(2) """ function dynamap(ps1,ps2,ps3,ps4) - #ps1 : time evolved system density matrix starting from initial state up - #ps2 : time evolved system density matrix starting from initial state down - #ps3 : time evolved system density matrix starting from initial state (up + down)/sqrt(2) - #ps4 : time evolved system density matrix starting from initial state (up - i*down)/sqrt(2) ϵ = [ ps1[1]-real(ps1[3])-imag(ps1[3]) ps2[1]-real(ps2[3])-imag(ps2[3]) ps3[1]-real(ps3[3])-imag(ps3[3]) ps4[1]-real(ps4[3])-imag(ps4[3]); ps1[4]-real(ps1[3])-imag(ps1[3]) ps2[4]-real(ps2[3])-imag(ps2[3]) ps3[4]-real(ps3[3])-imag(ps3[3]) ps4[4]-real(ps4[3])-imag(ps4[3]); @@ -453,6 +489,11 @@ function ttm2_evolve(numsteps, T, ps0) return reshape.(ps, 2, 2) end +""" + entropy(rho) + + +""" function entropy(rho) λ = eigen(rho).values return real(sum(map(x-> x==0 ? 0 : -x*log(x), λ))) @@ -469,6 +510,7 @@ end """ randisometry([T=Float64], dims...) + Construct a random isometry """ randisometry(T::Type, d1::Int, d2::Int) = d1 >= d2 ? Matrix(qr!(randn(T, d1, d2)).Q) : Matrix(lq!(randn(T, d1, d2)).Q) @@ -517,20 +559,3 @@ function MPOtoVector(mpo::MPO) H[N] = reshape(H[N], dims[1], 1, dims[2], dims[3]) return H end - -# Doesn't really work because indices change order in ITensors -function MPStoVector(mps::MPS) - N = length(mps) - A = [Array(mps[i], mps[i].inds...) for i=1:N] - dims = size(A[1]) - A[1] = reshape(A[1], 1, dims...) - A[1] = permutedims(A[1], [1,3,2]) - for i=2:N-1 - A[i] = permutedims(A[i], [3,1,2]) - end - dims = size(A[N]) - A[N] = reshape(A[N], dims..., 1) - A[N] = permutedims(A[N], [1,3,2]) - return A -end - diff --git a/src/models.jl b/src/models.jl index 3d3e0b6..329f792 100644 --- a/src/models.jl +++ b/src/models.jl @@ -522,24 +522,6 @@ function chaincoeffs_ohmic(nummodes, α, s; ωc=1, soft=false) end end -#Deprecated# -function getchaincoeffs(nummodes, α, s, beta, ωc=1) - matlabdir = ENV["MATDIR"] - astr = paramstring(α, 2) - bstr = paramstring(beta, 3) - datfname = string("jacerg_","a",astr,"s$s","beta",bstr,".csv") - chaincoeffs = readdlm(string(matlabdir,datfname),',',Float64,'\n') - es = ωc .* chaincoeffs[:,1] - ts = ωc .* chaincoeffs[1:end-1,2] - c0 = ωc * sqrt(chaincoeffs[end,2]/pi) - Nmax = length(es) - if Nmax < nummodes - throw(ErrorException("no data for nummodes > $Nmax")) - end - return [es[1:nummodes], ts[1:nummodes-1], c0] -end -## - """ readchaincoeffs(fdir, params...) @@ -714,11 +696,11 @@ Generate MPO for a pure dephasing model, defined by the Hamiltonian The spin is on site 1 of the MPS and the bath modes are to the right. ### Arguments - * `ΔE::Real`: energy splitting of the spin - * `dchain::Int`: physical dimension of the chain sites truncated Hilbert spaces - * `Nchain::Int`: number of sites in the chain - * `chainparams::Array{Real,1}`: chain parameters for the bath chain. The chain parameters are given in the standard form: `chainparams` ``=[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]``. - * `tree::Bool`: if true, return a `TreeNetwork` object, otherwise return a vector of MPO tensors +* `ΔE::Real`: energy splitting of the spin +* `dchain::Int`: physical dimension of the chain sites truncated Hilbert spaces +* `Nchain::Int`: number of sites in the chain +* `chainparams::Array{Real,1}`: chain parameters for the bath chain. The chain parameters are given in the standard form: `chainparams` ``=[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]``. +* `tree::Bool`: if true, return a `TreeNetwork` object, otherwise return a vector of MPO tensors """ function puredephasingmpo(ΔE, dchain, Nchain, chainparams; tree=false) u = unitmat(2) @@ -743,10 +725,10 @@ The interactions are nearest-neighbour, with the first N/2-1 sites corresponding # Arguments - * `N::Int`: number of sites in the chain - * `ϵd::Real`: energy of the impurity site at the center, as Ed - μ, where μ is the chemical potential - * chainparams1::Array{Real,1}: chain parameters for the first lead - * chainparams2::Array{Real,1}: chain parameters for the second lead +* `N::Int`: number of sites in the chain +* `ϵd::Real`: energy of the impurity site at the center, as Ed - μ, where μ is the chemical potential +* `chainparams1::Array{Real,1}`: chain parameters for the first lead +* `chainparams2::Array{Real,1}`: chain parameters for the second lead The chain parameters are given in the standard form: `chainparams` ``=[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]``. """ @@ -857,19 +839,19 @@ end """ interleaved_tightbinding_mpo(N, ϵd, chainparams1, chainparams2) - Generate MPO for a tight-binding chain of N fermionic sites with a single impurity site (fermionic as well) - of energy ϵd. The impurity is coupled to two leads, each described by a set of chain parameters. - The interactions are next-nearest-neighbour, with the first site corresponding to the impurity, and the - two chains organised in an interleaved fashion. +Generate MPO for a tight-binding chain of N fermionic sites with a single impurity site (fermionic as well) +of energy ϵd. The impurity is coupled to two leads, each described by a set of chain parameters. +The interactions are next-nearest-neighbour, with the first site corresponding to the impurity, and the +two chains organised in an interleaved fashion. - # Arguments +# Arguments - * `N::Int`: number of sites in the chain - * `ϵd::Real`: energy of the impurity site at the first site, as Ed - μ, where μ is the chemical potential - * chainparams1::Array{Real,1}: chain parameters for the first lead - * chainparams2::Array{Real,1}: chain parameters for the second lead +* `N::Int`: number of sites in the chain +* `ϵd::Real`: energy of the impurity site at the first site, as Ed - μ, where μ is the chemical potential +* chainparams1::Array{Real,1}: chain parameters for the first lead +* chainparams2::Array{Real,1}: chain parameters for the second lead - The chain parameters are given in the standard form: `chainparams` ``=[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]``. +The chain parameters are given in the standard form: `chainparams` ``=[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]``. """ function interleaved_tightbinding_mpo(N, ϵd, chainparams1, chainparams2) @@ -1333,31 +1315,32 @@ end """ protontransfermpo(ω0e,ω0k,x0e,x0k, Δ, dRC, d, N, chainparams, RCparams, λreorg) -Generate a MPO for a system described in space with a reaction coordinate (RC) tensor. The RC tensor is coupled to a bosonic bath, taking into account the induced reorganization energy. +Generate a MPO for a system described in space with a reaction coordinate (RC) tensor. The Hamiltonian of the two-level system and of the reaction coordinate tensor reads `` H_S + H_{RC} + H_{int}^{S-RC} = \\omega^0_{e} |e\\rangle \\langle e| + \\omega^0_{k} |k\\rangle \\langle k| + \\Delta (|e\\rangle \\langle k| + |k\\rangle \\langle e|) + \\omega_{RC} (d^{\\dagger}d + \\frac{1}{2}) + g_{e} |e\\rangle \\langle e|( d + d^{\\dagger})+ g_{k} |k \\rangle \\langle k|( d + d^{\\dagger}) `` +The RC tensor is coupled to a bosonic bath, taking into account the induced reorganization energy `` H_B + H_{int}^{RC-B} = \\int_{-∞}^{+∞} dk ω_k b_k^\\dagger b_k - (d + d^{\\dagger})\\int_0^∞ dω\\sqrt{J(ω)}(b_ω^\\dagger+b_ω) + \\lambda_{reorg}(d + d^{\\dagger})^2 -``. `` -\\lambda_{reorg} = \\int \\frac{J(\\omega)}{\\omega}d\\omega -``. - +with +`` +\\lambda_{reorg} = \\int \\frac{J(\\omega)}{\\omega}d\\omega. +`` # Arguments -* `ω0e`: enol energy at x=0 -* `ω0k`: keto energy at x=0 +* `ω0e`: enol energy at reaction coordinate value x=0 +* `ω0k`: keto energy at reaction coordinate value x=0 * `x0e`: enol equilibrium displacement * `x0k`: keto equilibrium displacement * `Δ`: direct coupling between enol and keto * `dRC`: fock space of the RC tensor * `d`: number of Fock states of the chain modes -* `N`: length of the chain +* `N`: length of the bosonic chain * `chainparams`: chain parameters, of the form `chainparams`=``[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]``, can be chosen to represent any arbitrary spectral density ``J(ω)`` at any temperature. -* `RCparams`: RC tensor parameter, of the form `RCparams`=``[ωRC,-g/x]`` +* `RCparams`: RC tensor parameter, of the form `RCparams`=``[ω_RC,-g/x]`` * `λreorg`: reorganization energy """ function protontransfermpo(ω0e,ω0k,x0e,x0k, Δ, dRC, d, N, chainparams, RCparams, λreorg) diff --git a/src/mpsBasics.jl b/src/mpsBasics.jl index 9c9187a..4d8bee6 100644 --- a/src/mpsBasics.jl +++ b/src/mpsBasics.jl @@ -1,10 +1,10 @@ """ - orthcentersmps(A) + orthcentersmps(A::Vector) Compute the orthoganality centres of MPS `A`. Return value is a list in which each element is the corresponding site tensor of `A` with the -orthoganility centre on that site. Assumes `A` is right normalised. +orthogonality centre on that site. Assumes `A` is right normalised. """ function orthcentersmps(A::Vector)