diff --git a/ChainOhmT/chainOhmT.jl b/ChainOhmT/chainOhmT.jl index 16f8ea3..cfcbc9c 100644 --- a/ChainOhmT/chainOhmT.jl +++ b/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 + diff --git a/ChainOhmT/fermionicT/chaincoeffs.h5 b/ChainOhmT/fermionicT/chaincoeffs.h5 new file mode 100644 index 0000000..b93f964 Binary files /dev/null and b/ChainOhmT/fermionicT/chaincoeffs.h5 differ diff --git a/ChainOhmT/mcdis2.jl b/ChainOhmT/mcdis2.jl index af80714..4ce8834 100644 --- a/ChainOhmT/mcdis2.jl +++ b/ChainOhmT/mcdis2.jl @@ -69,7 +69,7 @@ MCDIS Multiple-component discretization procedure. in multi-component form. """ -function mcdis(N,eps0,quad::Function,Nmax) +function mcdis(N,eps0,quad::Function,Nmax,idelta,mc,AB,wf,mp,irout) suc = true f = "Ncap exceeds Nmax in mcdis with irout = " @@ -104,7 +104,7 @@ function mcdis(N,eps0,quad::Function,Nmax) xwm = Array{Float64}(undef, mc*Ncap, 2) for i=1:mc im1tn = (i-1)*Ncap - xw = quad(Ncap,i,uv) + xw = quad(Ncap,i,uv,mc,AB,wf) xwm[im1tn+1:im1tn+Ncap,:] = xw end diff --git a/ChainOhmT/ohmicT/chaincoeffs.h5 b/ChainOhmT/ohmicT/chaincoeffs.h5 new file mode 100644 index 0000000..0ce402f Binary files /dev/null and b/ChainOhmT/ohmicT/chaincoeffs.h5 differ diff --git a/ChainOhmT/ohmicT/chaincoeffs_ohmic_a0.01wc0.035xc0.035beta100.csv b/ChainOhmT/ohmicT/chaincoeffs_ohmic_a0.01wc0.035xc0.035beta100.csv new file mode 100644 index 0000000..56599c0 --- /dev/null +++ b/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 diff --git a/ChainOhmT/quadohmT.jl b/ChainOhmT/quadohmT.jl index 9851c4f..6afed43 100644 --- a/ChainOhmT/quadohmT.jl +++ b/ChainOhmT/quadohmT.jl @@ -1,73 +1,96 @@ -function quadfinT(N,i,uv) +function quadfinT(N,i,uv,mc,AB,wf) if (i>1 && ix!=0, xw[:,2]) #index = min(find(xw(:,2)~=0)) xw = xw[index:length(xw[:,2]),:] - Ibis = sortper(xw[:,1]) #xw = sortrows(xw,1) + Ibis = sortperm(xw[:,1]) #xw = sortrows(xw,1) xw = xw[Ibis,:] Ncap = size(xw,1) @@ -27,7 +27,8 @@ function stieltjes(N,xw) #return ab error("N in sti out of range") end - s0 = ones(1,Ncap)*xw[:,2] + s0 = (ones(1,Ncap)*xw[:,2])[1] + ab = zeros(N, 2) ab[1,1] = transpose(xw[:,1])*xw[:,2]/s0 ab[1,2] = s0 if N==1 @@ -50,23 +51,22 @@ function stieltjes(N,xw) #return ab # % end # % - c = 1e240 - p2 = c*p2 - s0 = c^2*s0 + #c = 1e240 + #p2 = c*p2 + #s0 = c^2*s0 for k=1:N-1 p0 = p1 p1 = p2 - p2 = (xw[:,1] - ab[k,1]).*p1 - ab[k,2]*p0 - s1 = transpose(xw[:,2])*(p2.^2) - s2 = transpose(xw[:,1])*(xw[:,2].*(p2.^2)) - - if (max(abs(p2))>huge)||(abs(s2)>huge) - error("impending overflow in stieltjes for k=$(k)") - end - if abs(s1)huge)||(abs(s2)>huge) + # error("impending overflow in stieltjes for k=$(k)") + # end + # if abs(s1)", "Thibaut Lacroix "] -version = "1.0" +version = "1.0.0" [deps] Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" @@ -10,17 +10,19 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JLD = "4138dd39-2aa7-5051-a626-17a0bb65d9c8" +Jacobi = "83f21c0b-4282-5fbc-9e3f-f6da3d2e584c" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Preferences = "21216c6a-2e73-6563-6e65-726566657250" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" -QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" -Jacobi = "83f21c0b-4282-5fbc-9e3f-f6da3d2e584c" [compat] TensorOperations = "4.0.7" diff --git a/README.md b/README.md index 55c1fb6..c5b5119 100644 --- a/README.md +++ b/README.md @@ -1,20 +1,23 @@ +MPSDynamics.jl logo + # MPSDynamics.jl *Tensor network simulations for finite temperature, open quantum system dynamics.* [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.5106435.svg)](https://doi.org/10.5281/zenodo.5106435) [![license](https://img.shields.io/badge/License-GPL_3.0-orange.svg)](https://github.com/angusdunnett/MPSDynamics/blob/master/LICENSE) [![documentation workflow](https://github.com/angusdunnett/MPSDynamics/actions/workflows/docs.yml/badge.svg)](https://angusdunnett.github.io/MPSDynamics/) + This package is intended to provide an easy to use interface for performing tensor network simulations on Matrix Product States (MPS). MPSDynamics.jl is a versatile package which supports both chain and (loop-free) tree MPS, as well as providing a choice of several time evolution algorithms. The package also provides strong support for the measurement of observables, as well as the storing and logging of data, which makes it a useful tool for the study of many-body physics. The package was originally developed with the aim of studying open system dynamics at finite temperature using -the T-TEDOPA mapping [1], however the methods implemented can equally be applied to other areas of physics. +the T-TEDOPA mapping [^1], however the methods implemented can equally be applied to other areas of physics. The methods currently implemented are -* 1-site TDVP on tree and chain MPS [2] -* 2-site TDVP on chain MPS [2] -* a variant of 1-site TDVP with dynamic bond-dimensions on chain MPS [3] +* 1-site TDVP on tree and chain MPS [^2] +* 2-site TDVP on chain MPS [^2] +* a variant of 1-site TDVP with dynamic bond-dimensions on chain MPS [^3] The elementary tensor operations are implemented in all cases using the [TensorOperations.jl](https://github.com/Jutho/TensorOperations.jl) package. @@ -149,6 +152,8 @@ Import["~/MPSDynamics/XXXXX/dat_XXXXX.jld",{"HDF5","Datasets","/data/sz"}] # Publications Publications which make use of MPSDynamics: +* Lacroix et al. From Non-Markovian Dissipation to Spatiotemporal Control of Quantum Nanodevices. *Quantum* 8, 1305, April 2024 + * [https://doi.org/10.22331/q-2024-04-03-1305](https://doi.org/10.22331/q-2024-04-03-1305) * Riva et al. Thermal cycle and polaron formation in structured bosonic environments. *Phys. Rev. B* 108(19):195138, November 2023 * [https://doi.org/10.1103/PhysRevB.108.195138](https://doi.org/10.1103/PhysRevB.108.195138) @@ -186,9 +191,9 @@ Contributions are welcome! Don't hesitate to contact us if you # References -* [[1]](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.123.090402) D. Tamascelli, A. Smirne, J. Lim, S. F. Huegla, and M. B. Plenio, Physical Review Letters 123, 090402 (2019) arXiv: 1811.12418 +[^1]: [D. Tamascelli, A. Smirne, J. Lim, S. F. Huegla, and M. B. Plenio, Physical Review Letters 123, 090402 (2019) arXiv: 1811.12418](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.123.090402) -* [[2]](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.94.165116) J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Physical Review B 94, 165116 (2016), arXiv: 1408.5056 +[^2]: [J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken, and F. Verstraete, Physical Review B 94, 165116 (2016), arXiv: 1408.5056](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.94.165116) -* [[3]](https://doi.org/10.1103/PhysRevB.104.214302) A. J. Dunnett & A. W. Chin, Physical Review B, 104(21), 214302 (2021) +[^3]: [A. J. Dunnett & A. W. Chin, Physical Review B, 104(21), 214302 (2021)](https://doi.org/10.1103/PhysRevB.104.214302) diff --git a/examples/anderson_model_double.jl b/examples/anderson_model_double.jl new file mode 100644 index 0000000..91028ac --- /dev/null +++ b/examples/anderson_model_double.jl @@ -0,0 +1,150 @@ +using MPSDynamics, Plots, LaTeXStrings +import MPSDynamics: productstatemps, tightbinding_mpo, mpsembed! + +#---------------------------- +# Physical parameters +#---------------------------- + +N = 40 # number of chain sites +β = 10.0 # inverse temperature +μ = 0. # chemical potential +Ed = 0.3 # energy of the impurity +ϵd = Ed - μ # energy of the impurity minus the chemical potential + +#----------------------------------------- +# Dispersion relation and chain parameters +#----------------------------------------- + +function ϵ(x) + return x +end + +function J(x) + return sqrt(1 - x^2) # semi-circular density of states +end + +chainparams1 = chaincoeffs_fermionic(N, β, 1.0; ϵ, J, save=false) # empty +chainparams2 = chaincoeffs_fermionic(N, β, 2.0; ϵ, J, save=false) # filled + +#= +curdir = @__DIR__ +dir_chaincoeff = abspath(joinpath(curdir, "../ChainOhmT/fermionicT/")) +chainparams2 = readchaincoeffs("$dir_chaincoeff/chaincoeffs.h5", β, 2.0) #filled +chainparams1 = readchaincoeffs("$dir_chaincoeff/chaincoeffs.h5", β, 1.0) #empty +=# + +#----------------------- +# Simulation parameters +#----------------------- + +dt = 0.25 # time step +T = 15.0 # simulation time +method = :DTDVP # time-evolution method +Dmax = 150 # MPS max bond dimension +prec = 0.00001 # precision for the adaptive TDVP + +#--------------------------- +# MPO and initial state MPS +#--------------------------- + +H = tightbinding_mpo(N, ϵd, chainparams1, chainparams2) + +ψ = unitcol(2,2) # (0,1) filled impurity state +A = productstatemps(physdims(H), state=[fill(unitcol(2,2), N)..., ψ, fill(unitcol(1,2), N)...]) # MPS +mpsembed!(A, 2) + +#---------------------------------------------------- +# Definition of observables for the double chain +#---------------------------------------------------- + +# For the double chain +ob1 = OneSiteObservable("chain1_filled_occup", numb(2), (1,N)) +ob2 = OneSiteObservable("chain2_empty_occup", numb(2), (N+2, 2N+1)) +ob3 = OneSiteObservable("system_occup", numb(2), N+1) + +#------------- +# Simulation +#------------ + +A, dat = runsim(dt, T, A, H; + name = "Anderson impurity problem (folded chain)", + method = method, + obs = [ob1, ob2, ob3], + convobs = [ob1], + params = @LogParams(N, ϵd, β, c1, c2), + convparams = [prec], + Dlim = Dmax, + savebonddims = true, # we want to save the bond dimension + verbose = false, + save = false, + plot = true, + ); + +#---------------- +# Post-processing +#---------------- + +# Reshaping the vector to a column matrix and horizontal concatenation +system_occup_col = reshape(dat["data/system_occup"], :, 1) +occ = hcat(dat["data/chain1_filled_occup"]', system_occup_col) +occ = vcat(occ', dat["data/chain2_empty_occup"]) + +#------------- +# Plots +#------------- + +# Plot the system occupation +p1 = plot( + dat["data/times"], + dat["data/system_occup"], + xlabel = L"$t$", + ylabel = L"$n_d$", + title = "System Occupation", + size = (700, 500) +) + +# Plot the occupation of the chain sites +p2 = heatmap( + collect(1:2*N+1), + dat["data/times"], + transpose(occ), # Use the matrix form + cmap = :coolwarm, + aspect_ratio = :auto, + xlabel = L"$N_{i,j}$ chain sites", + ylabel = L"$t$", + title = "Chain Occupation", + colorbar = true, + size = (700, 500) +) + +# Plot the bond dimensions +p3 = heatmap( + collect(1:2*N+2), + dat["data/times"], + transpose(dat["data/bonddims"]), + cmap = :magma, + aspect_ratio = :auto, + xlabel = L"$N_{i,j}$ chain sites", + ylabel = L"$t$", + title = "Bond Dimensions", + colorbar = true, + size = (700, 500) +) + +# Define indices for columns to be plotted +columns_to_plot = [1, 5, 10, 15, 20] + +# Plot vertical slices for occupancy +p4 = plot(title = "Chain occupation") +for col in columns_to_plot + plot!(p4, occ[:, col], label = L"$t =$"*"$col", xlabel = L"$N_{i,j}$ chain sites", ylabel = "chain occupation") +end + +# Plot vertical slices for bond dimensions +p5 = plot(title = "Bond Dimensions") +for col in columns_to_plot + plot!(p5, dat["data/bonddims"][:, col], label = L"$t =$"*"$col", xlabel = L"$N_{i,j}$ chain sites", ylabel = L"$\chi$") +end + +# Display the plots +plot(p2, p3, p4, p5, p1, layout = (3, 2), size = (1400, 1200)) diff --git a/examples/anderson_model_interleaved.jl b/examples/anderson_model_interleaved.jl new file mode 100644 index 0000000..e845cd5 --- /dev/null +++ b/examples/anderson_model_interleaved.jl @@ -0,0 +1,180 @@ +using MPSDynamics, Plots, LaTeXStrings +import MPSDynamics: mpsembed!, interleaved_tightbinding_mpo, productstatemps + +#---------------------------- +# Physical parameters +#---------------------------- + +N = 40 # number of chain sites +β = 10.0 # inverse temperature +μ = 0. # chemical potential +Ed = 0.3 # energy of the impurity +ϵd = Ed - μ # energy of the impurity minus the chemical potential + +#----------------------------------------- +# Dispersion relation and chain parameters +#----------------------------------------- + +function ϵ(x) + return x +end + +function J(x) + return sqrt(1 - x^2) # semi-circular density of states +end + +chainparams1 = chaincoeffs_fermionic(N, β, 1.0; ϵ, J, save=false) # empty +chainparams2 = chaincoeffs_fermionic(N, β, 2.0; ϵ, J, save=false) # filled + +#= +curdir = @__DIR__ +dir_chaincoeff = abspath(joinpath(curdir, "../ChainOhmT/fermionicT/")) +chainparams2 = readchaincoeffs("$dir_chaincoeff/chaincoeffs.h5", β, 2.0) #filled +chainparams1 = readchaincoeffs("$dir_chaincoeff/chaincoeffs.h5", β, 1.0) #empty +=# + +#----------------------- +# Simulation parameters +#----------------------- + +dt = 0.25 # time step +T = 15.0 # simulation time +method = :DTDVP # time-evolution method +Dmax = 150 # MPS max bond dimension +prec = 0.0001 # precision for the adaptive TDVP + +#--------------------------- +# MPO and initial state MPS +#--------------------------- + +H = interleaved_tightbinding_mpo(N, ϵd, chainparams1, chainparams2) + +ψ = unitcol(2,2) # (0,1) filled impurity state +Tot = Any[] +push!(Tot, ψ) +for i in 1:(N) + push!(Tot, unitcol(2,2)) + push!(Tot, unitcol(1,2)) +end + +A = productstatemps(physdims(H), state=Tot) # MPS +mpsembed!(A, 2) # embed the MPS in a manifold of bond dimension 2 + + +#---------------------------------------------------- +# Definition of observables for the interleaved chain +#---------------------------------------------------- + +ob1 = OneSiteObservable("system_occup", numb(2), 1) +ob2 = OneSiteObservable("folded_chain_occup", numb(2), (2,2N+1)) + +#------------- +# Simulation +#------------ + +A, dat = runsim(dt, T, A, H; + name = "Anderson impurity problem (folded chain)", + method = method, + obs = [ob1, ob2], + convobs = [ob1], + params = @LogParams(N, ϵd, β, c1, c2), + convparams = [prec], + Dlim = Dmax, + savebonddims = true, # we want to save the bond dimension + verbose = false, + save = false, + plot = true, + ); + + +#------------------------------------------------------------- +# Post-processing of the data: unfolding the chain for clarity +#------------------------------------------------------------- +unfolded_occ = Vector{Vector{Float64}}() # Assuming the elements are of type Float64 +unfolded_bonds = Vector{Vector{Float64}}() # Adjust the type based on actual data + +# Populate unfolded_occ by iterating in the specific order mentioned +for i in 1:N # Adjusted for 1-based indexing + push!(unfolded_occ, dat[ "data/folded_chain_occup"][2N + 1 - 2i, :]) +end + +push!(unfolded_occ, dat["data/folded_chain_occup"][1,:]) + +for i in 2:N + push!(unfolded_occ, dat["data/folded_chain_occup"][2i,:]) +end + +# Populate unfolded_bonds similarly +for i in 1:(N+1) # Adjusted for 1-based indexing + push!(unfolded_bonds, dat["data/bonddims"][2N + 3 - 2i,:]) # Assuming bonddims is directly accessible +end + +push!(unfolded_bonds, dat["data/bonddims"][1,:]) + +for i in 2:(N+1) + push!(unfolded_bonds, dat["data/bonddims"][2i,:]) +end + +unfolded_bonds_matrix = hcat(unfolded_bonds...)' +unfolded_occ_matrix = hcat(unfolded_occ...)' + +#------------- +# Plots +#------------- + +# Plot the system occupation +p1 = plot( + dat["data/times"], + dat["data/system_occup"], + xlabel = L"$t$", + ylabel = L"$n_d$", + title = "System Occupation", + size = (700, 500) +) + +# Plot the occupation of the chain sites +p2 = heatmap( + collect(1:2*N), + dat["data/times"], + transpose(unfolded_occ_matrix), # Use the matrix form + cmap = :coolwarm, + aspect_ratio = :auto, + xlabel = L"$N_{i,j}$ chain sites", + ylabel = L"$t$", + title = "Chain Occupation", + colorbar = true, + size = (700, 500) +) + +# Plot the bond dimensions +p3 = heatmap( + collect(1:2*N+2), + dat["data/times"], + transpose(unfolded_bonds_matrix), + cmap = :magma, + aspect_ratio = :auto, + xlabel = L"$N_{i,j}$ chain sites", + ylabel = L"$t$", + title = "Bond Dimensions", + colorbar = true, + size = (700, 500) +) + +# Define indices for columns to be plotted +columns_to_plot = [1, 5, 10, 15, 20] + + +# Plot vertical slices for occupancy +p4 = plot(title = "Chain occupation") +for col in columns_to_plot + plot!(p4, unfolded_occ_matrix[:, col], label = L"$t =$"*"$col", xlabel = L"$N_{i,j}$ chain sites", ylabel = "chain occupation") +end + +# Plot vertical slices for bond dimensions +p5 = plot(title = "Bond Dimensions") +for col in columns_to_plot + plot!(p5, unfolded_bonds_matrix[:, col], label = L"$t =$"*"$col", xlabel = L"$N_{i,j}$ chain sites", ylabel = L"$\chi$") +end + +# Display the plots +plot(p2, p3, p4, p5, p1, layout = (3, 2), size = (1400, 1200)) \ No newline at end of file diff --git a/examples/bath-observables.jl b/examples/bath-observables.jl new file mode 100644 index 0000000..caf0f39 --- /dev/null +++ b/examples/bath-observables.jl @@ -0,0 +1,173 @@ +using MPSDynamics, Plots, LaTeXStrings, QuadGK, LinearAlgebra, Interpolations, Revise + +import MPSDynamics: measuremodes, measurecorrs, mpsembed!, eigenchain, physical_occup + + +const ∞ = Inf +#---------------------------- +# Physical parameters +#---------------------------- + +d = 10 # number of Fock states of the chain modes +N = 60 # length of the chain +α = 0.01 # coupling strength +ω0 = 0.2 # TLS gap +s = 1 # ohmicity +ωc = 1. # Cut-off of the spectral density J(ω) +β = 20 # Thermalized environment +#β = ∞ # Case zero temperature T=0, β → ∞ + + +#----------------------- +# Simulation parameters +#----------------------- + +method = :TDVP1 # time-evolution method +conv = 3 # bond dimension for the TDVP1 +dt = 0.5 # time step +tfinal = 50.0 # simulation time + +#---------------------------- +# Ohmic spectral density +#---------------------------- + +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 + +#--------------------------- +# MPO and initial state MPS +#--------------------------- + +H = puredephasingmpo(ω0, 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) +ob2 = OneSiteObservable("sx", sx, 1) +ob3 = OneSiteObservable("chain_mode_occupation", numb(d), (2,N+1)) +ob4 = OneSiteObservable("c", crea(d), collect(2:N+1)) +ob5 = OneSiteObservable("cdag", crea(d), collect(2:N+1)) +ob6 = TwoSiteObservable("cdagc", crea(d), anih(d), collect(2:N+1), collect(2:N+1)) +ob7 = TwoSiteObservable("cdagcdag", crea(d), crea(d), collect(2:N+1), collect(2:N+1)) +ob8 = TwoSiteObservable("cc", anih(d), anih(d), collect(2:N+1), collect(2:N+1)) + +#------------- +# Simulation +#------------ + +A, dat = runsim(dt, tfinal, A, H, prec=1E-4; + name = "Bath observables in the pure dephasing model", + method = method, + obs = [ob1, ob2, ob3, ob4, ob5, ob6, ob7, ob8], + convobs = [ob1], + params = @LogParams(ω0, N, d, α, s, ψ), + convparams = conv, + reduceddensity = true, + verbose = false, + save = false, + plot = true, + ); + +#--------------------------- +# Post-processing +#--------------------------- + +T = length(dat["data/times"]) + +constr = Array{ComplexF64}(undef, N, N, T) +destr = Array{ComplexF64}(undef, N, N, T) +for t in 1:T + constr[:,:,t] = diagm(0 => dat["data/cdag"][:,t]) + destr[:,:,t] = diagm(0 => dat["data/c"][:,t]) +end + +omeg = eigenchain(cpars, nummodes=N).values +bath_occup = mapslices(X -> measuremodes(X, cpars[1], cpars[2]), dat["data/cdagc"], dims = [1,2]) +cdag_average = mapslices(X -> measuremodes(X, cpars[1], cpars[2]), constr, dims = [1,2]) +c_average = mapslices(X -> measuremodes(X, cpars[1], cpars[2]), destr, dims = [1,2]) +cc_average = mapslices(X -> measurecorrs(X, cpars[1], cpars[2]), dat["data/cc"], dims = [1,2]) +cdagcdag_average = mapslices(X -> measurecorrs(X, cpars[1], cpars[2]), dat["data/cdagcdag"], dims = [1,2]) + +correlations_c = [ + cc_average[i, j, t] - c_average[i, 1, t] .* c_average[j, 1, t] + for i in 1:size(cc_average, 1), j in 1:size(cc_average, 2), t in 1:size(cc_average, 3) +] +correlations_cdag = [ + cdagcdag_average[i, j, t] - cdag_average[i, 1, t] .* cdag_average[j, 1, t] + for i in 1:size(cdagcdag_average, 1), j in 1:size(cdagcdag_average, 2), t in 1:size(cdagcdag_average,3) +] + +bath_occup_phys = physical_occup(correlations_cdag[:,:,T], correlations_c[:,:,T], omeg, bath_occup[:,:,T], β, N) + +#-------------------- +# Analytical results +#-------------------- + +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)) + + +α_theo = 0.25 * α +function Jtherm(x) + if 1 >= x >= 0 + return +α_theo * abs(x)^s * (1 + coth(β*0.5*x)) + elseif -1 <= x <= 0 + return -α_theo * abs(x)^s * (1 + coth(β*0.5*x)) + else + return 0 + end +end + +bath_occup_analytical(ω, t) = abs(Jtherm(ω))/(ω^2)*2*(1-cos(ω*t)) + +#------------- +# Plots +#------------ + +ρ12 = abs.(dat["data/Reduced ρ"][1,2,:]) + +p1 = 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) +p1 = plot!(dat["data/times"], ρ12, lw=4, ls=:dash, label="Numerics") + +cumul = [bath_occup_analytical(omeg[i], tfinal)*(omeg[i+1]-omeg[i]) for i in 1:(length(omeg)-1)] + +p2 = plot(omeg[1:length(omeg)-1], cumul, lw = 4, linecolor=:black, + xlabel=L"\omega", ylabel=L"\langle n^b_\omega \rangle", label="Analytics", + title="Mode occupation in the extended bath") +p2 = plot!(omeg, bath_occup[:, :, T], lw=4, ls=:dash, label="Numerics") + +p3 = heatmap(omeg, omeg, abs.(real.(correlations_cdag[:,:,T]) .+ im*imag.(correlations_cdag[:,:,T])), + xlabel=L"\omega", + ylabel=L"\omega", title="Environmental correlations") + + +Mhalf = Int(length(omeg)*0.5)+1 +M = length(omeg) + +p4 = plot(omeg[Mhalf:M], bath_occup_phys, lw=4, + xlabel=L"\omega", ylabel=L"\langle n^b_\omega \rangle", + title="Mode occupation in the physical bath") + +plot(p1, p2, p3, p4, layout = (2, 2), size = (1400, 1200)) diff --git a/examples/plots/chain_occ_bet=2_double.pdf b/examples/plots/chain_occ_bet=2_double.pdf new file mode 100644 index 0000000..768b0f2 Binary files /dev/null and b/examples/plots/chain_occ_bet=2_double.pdf differ diff --git a/examples/plots/chain_occ_bet=2_folded.pdf b/examples/plots/chain_occ_bet=2_folded.pdf new file mode 100644 index 0000000..64fe576 Binary files /dev/null and b/examples/plots/chain_occ_bet=2_folded.pdf differ diff --git a/examples/plots/system_folded_beta2.pdf b/examples/plots/system_folded_beta2.pdf new file mode 100644 index 0000000..00b6831 Binary files /dev/null and b/examples/plots/system_folded_beta2.pdf differ diff --git a/examples/protontransfer.jl b/examples/protontransfer.jl new file mode 100644 index 0000000..1a85c2b --- /dev/null +++ b/examples/protontransfer.jl @@ -0,0 +1,273 @@ +#= + Example of a Proton Transfer model at zero temperature for an isolated system or with an environment characterised by 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. + + The RC tensor is initially displaced at the value γ corresponding to a space coordinate. + + The initial electronic system is initialized in the adiabatic LOW surface at the RC displacement. + + H_S + H_RC + H_int^{S-RC} = ω^0_{e} |e⟩ ⟨e| + ω^0_{k} |k⟩ ⟨k| + \\Delta (|e⟩ ⟨k| + |k⟩ ⟨ e|) + ω_{RC} (d^{\dagger}d + \\frac{1}{2}) + g_{e} |e⟩ ⟨e|( d + d^{\dagger})+ g_{k} |k⟩ ⟨k|( d + d^{\dagger}) +`` + H_B + H_int^{RC-B} = ∫_{-∞}^{+∞} dk ω_k b_k^\dagger b_k - (d + d^{\dagger})∫_0^∞ dω\\sqrt{J(ω)}(b_ω^\\dagger+b_ω) + λ_{reorg}(d + d^{\\dagger})^2 +``. + λ_{reorg} = ∫ \frac{J(ω)}{ω}dω + +=# + +using MPSDynamics, Plots, LaTeXStrings, ColorSchemes, PolyChaos, LinearAlgebra + +#---------------------------- +# Physical parameters +#---------------------------- +# Enol / Keto + +ω0e= 0.8 # Enol energy at x=0 +ω0k= 0.8 # Keto energy at x=0 + +x0e = -0.25 # Enol Equilibrium displacement +x0k = 0.25 # Keto Equilibrium displacement + +Δ = 0.05 # Direct coupling between enol and keto + +m=1.83E3 # Mass of the reaction coordinate particle. 1.83E3 is the proton mass in atomic units. +ħ=1 # Atomic Units convention + +dFockRC = 25 # Fock space of the RC tensor + +ωRC = 0.0347 # Frequency of the RC tensor + +γ = sqrt(m*ωRC/2)*x0e # γ = displacement RC ; γ = \sqrt{m\omega/2} x_disp + +cparsRC = [ωRC,ωRC*sqrt(m*ωRC/2)] # Energy and g RC coupling parameter + +isolated = true # isolated = true : no environment ; isolated = false : system coupled to a bosonic environment +# Creates the chain depending on the isolated condition. The parameters for isolated = false can be modified as desired. +if isolated + d=1; N=2; α = 0.0 # number of Fock states of the chain modes ; length of the chain ; coupling strength +else + d=4; N=40; α = 0.008 # number of Fock states of the chain modes ; length of the chain ; coupling strength +end + +s = 1 # ohmicity + +ωc = 2*ωRC # Cut-off of the spectral density J(ω) + +λreorg = (2*α*ωc)/s # Reorganisation Energy taken into account in the Hamiltonian + +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 = 10.0 # time step + +tfinal = 2000.0 # simulation time + +numsteps = length(collect(0:dt:tfinal))-1 + +method = :TDVP1 # time-evolution method + +D = 5 # MPS bond dimension + +#----------------------- +# Plot parameters +#----------------------- + +xlist_rho =collect(-1.2:0.05:1.2) # x definition for the reduced density matrix expressed in space +# other values of xlist_rho can be chosen to gain numerical time + +palette = ColorSchemes.okabe_ito # color palette for plots + +#--------------------------- +# MPO and initial state MPS +#--------------------------- + +# Initial electronic system in the left well, at mean displacement . It results to a superposition of 1 and 2 +# This state is obtained with the diagonalization of the diabatic Hamiltonian at x=γ +ψ = zeros(2) +X2 = zeros(2) + +Ae = ω0e - 0.5*m*ωRC^2*(x0e)^2 +Ak = ω0k - 0.5*m*ωRC^2*(x0k)^2 +xini = γ/(sqrt(m*ωRC/2)) + +Adumb = (ω0e-ω0k+0.5*m*ωRC^2*((xini-x0e)^2-(xini-x0k)^2)) +Bdumb = sqrt((Adumb)^2 + 4*Δ^2) +X2[1] = (Adumb-Bdumb)/(2*Δ) +X2[2] = 1 + +mod_X2 = sqrt(sum(X2[i]^2 for i=1:2)) + +ψ[1] = X2[1] / mod_X2 +ψ[2] = X2[2] / mod_X2 + +#### Coherent state . Initialise the RC oscillator at displacement γ thanks to the Taylor expression of the displacement operator + +coherent_RC = [fill(unitcol(1,dFockRC), 1)...] +chainpop_RCini =zeros(dFockRC,2) + +chainpop_RCini[1,2]=1 +chainpop_RCini[1,1]=exp(-(γ^2)/2) + +for j=2:dFockRC + chainpop_RCini[j,1]=chainpop_RCini[j-1,1]*(γ/sqrt(j-1)) +end + +coherent_RC[1] = hcat(chainpop_RCini[:,1]...) #Manipulation to get good format for MPS, does not change value + +H = protontransfermpo(ω0e, ω0k, x0e, x0k, Δ, dFockRC, d, N, cpars, cparsRC, λreorg) + +A = productstatemps(physdims(H), state=[ψ, coherent_RC..., fill(unitcol(1,d), N)...]) # MPS representation of |ψ>|Displaced RC>|Vacuum> + +Eini = measurempo(A,H) +print("\n Energy =",Eini) + +#--------------------------- +# Definition of observables +#--------------------------- + +ob1 = OneSiteObservable("sz", sz, 1) +#------------- +ob2 = OneSiteObservable("RC displacement", MPSDynamics.disp(dFockRC,ωRC,m), 2) + +#------------- +# Simulation +#------------ + +A, dat = runsim(dt, tfinal, A, H; + name = "proton transfer and tunneling at zero temperature", + method = method, + obs = [ob1, ob2], + convobs = [ob1, ob2], + params = @LogParams(N, d, α, s), + convparams = D, + Nrho = 2, # Need to precise that the reduced density matrix is + # calculated for the 2nd tensor. The value by default is 1. + reduceddensity= true, + verbose = false, + save = true, + plot = false, + ); + + +#------------- +# Plots +#------------ + +#### Plot the doublewell in the adiabatic basis. It is obtained with the diagonalization of the diabatic space Hamiltonian ##### + +function fp(ɸ) + H11 = Ae + 0.5*m*ωRC^2*(ɸ-x0e)^2 + H22 = Ak + 0.5*m*ωRC^2*(ɸ-x0k)^2 + H12 = Δ + λp = ((H11+H22)/2)+0.5*sqrt((H11-H22)^2+4*H12^2) + return λp +end + +function fm(ɸ) + H11 = Ae + 0.5*m*ωRC^2*(ɸ-x0e)^2 + H22 = Ak + 0.5*m*ωRC^2*(ɸ-x0k)^2 + H12 = Δ + λm = ((H11+H22)/2)-0.5*sqrt((H11-H22)^2+4*H12^2) + return λm +end + +function gaussian(x) + return Eini .+((1/0.06*sqrt(2*pi)).*exp.(-(x.+(0.25)).^2 ./(2*0.06^2))) .*0.001 +end + +ɸ = [i for i in -0.75:0.001:0.75] +ɸwp = [i for i in -0.45:0.001:-0.05] +xmin = -0.75 ; xmax = 0.75 ; ymin = 0.7 ; ymax = 0.9 + +plt=plot(ɸ,fm, label="LOW",xlabel="Reaction Coordinate (arb. units)",ylabel=L"$\omega$ (arb. units)",linewidth=4,thickness_scaling = 1, bg_legend=:transparent, legendfont=16, legendfontsize=16, xguidefontsize=16, yguidefontsize=16, tickfontsize=(12), xlims=(xmin,xmax),ylims=(ymin,ymax),xticks=(xmin:0.25:xmax),yticks=(ymin:0.05:ymax),legend=(0.6,0.6), grid=false,linecolor =palette[8]) + +plt=plot!(ɸ,fp,label="UP",lc =palette[3], linewidth=3,thickness_scaling = 1) + +plt = hline!([Eini], lw=2, lc=:black, ls=:dash, label=L"\omega_{system}") + +plt=plot!(ɸwp,gaussian,linecolor=:black,linewidth=3, label="", fillrange = (Eini , gaussian(ɸwp)),c = :black) + +display(plt) + +##### Results ##### + +display(plot(dat["data/times"], dat["data/RC displacement"],label=" (arb. units)", linecolor =:black, xlabel="Time (arb. units)",ylabel="", title="", linewidth=4, legend=:none, tickfontsize=10,xlims=(0,2000))) + +#### Reduced Density Matrix in space##### + +# Loads the Hermite coefficients to calculate the eigenstates of the harmonic oscillator. The coefficient 2 comes from the definition of the physics field. +αherm, βherm = 2.0 .*rm_hermite(dFockRC) + +#Fct building the nth harmonic eigenstate at coordinate x +# sqrt(m*ωRC/ħ)*sqrt(2)*x in the Hermite function to recover the form of the physical Hermite polynomial in the atomic units +function Ψ(n,x) +return Acoef(n)*(m*ωRC/(ħ*π))^(1/4)*exp(-m*ωRC*(x)^2/(2*ħ))*PolyChaos.evaluate(n,(sqrt(m*ωRC/ħ)*sqrt(2)*x),αherm,βherm) +end + + +#Function to compute the pre factor of Ψ(n,x) + function Acoef(n) + if n == 0 + return 1.0 + else + return (1/sqrt(n*2))*2^(1/2)*Acoef(n-1) + end +end + +#### Calculation of functions to express the reduced densty matrix in space. These functions build the matrix of oscillator eigenstates +# Ψ(i,x) Ψ(j,x) ∀i,j ∈ dFockRC^2 and ∀ x ∈ xlist_rho +Pmatrix = zeros(length(xlist_rho),dFockRC) +Pmatrixnormed = zeros(length(xlist_rho),dFockRC) +normalization = zeros(dFockRC) + +for n=1:dFockRC + for x=1:length(xlist_rho) + Pmatrix[x,n] = Ψ(n-1,xlist_rho[x]) + end +end + +for n=1:dFockRC + normalization[n] = sqrt(sum(Pmatrix[k,n]^2 for k=1:length(xlist_rho))) +end + + +for n=1:dFockRC + Pmatrixnormed[:,n] = Pmatrix[:,n]/normalization[n] +end + +### Converts the reduced density matrix expressed in Fock states in a reduced density matrix expressed in space. +ρx=zeros(ComplexF64,length(xlist_rho),length(xlist_rho),length(dat["data/times"])) +for t=1:length(dat["data/times"]) + print(" \n t = ", dat["data/times"][t]) + for x=1:length(xlist_rho) + for y=1:length(xlist_rho) + for n=1:dFockRC + for m=1:dFockRC + ρx[x,y,t] += Pmatrixnormed[x,n]*Pmatrixnormed[y,m]*dat["data/Reduced ρ"][n,m,t] + end + end + end + end +end + +#Creates a GIF of the diagonal elements of the reduced density over time, ie the population dynamics. Here, the reduced density matrix is expressed in space. + +anim = @animate for t=1:length(dat["data/times"]) + k=(dat["data/times"][t]) + print("\n k = ", k) + plot(xlist_rho, abs.(diag(ρx[:,:,t])), title="Time = $k (arb. units)", xlabel=L"Diag( $||\rho_{RC}(x,x)||$ )", ylabel="Amplitude", xlims=(-0.5, 0.5), ylims=(0,0.25),legend=:none) +end +display(gif(anim, "gif_population.gif", fps = 2.5)) + +#Creates a GIF of the reduced density expressed over time. Here, the reduced density matrix is expressed in space. Diagonal elements are populations +#whereas anti-diagonal elements represent coherences + +anim = @animate for t=1:length(dat["data/times"]) + k=(dat["data/times"][t]) + print("\n k = ", k) + (plot(heatmap(xlist_rho, xlist_rho, abs.(ρx[:,:,t]), c=:Blues ), title="Time = $k (arb. units)", xlabel=L"$x$ (arb. units)",ylabel=L"$x$ (arb. units)", tickfontsize=(12),colorbar_title = L"||\rho_{RC}(x,x)||", legend=:none, xlims=(-0.5, 0.5), ylims=(-0.5,0.5), clims=(0,0.18),aspect_ratio=:equal)) +end +display(gif(anim, "gif_reducedrho.gif", fps = 2.5)) diff --git a/examples/puredephasing.jl b/examples/puredephasing.jl new file mode 100644 index 0000000..ce8c0db --- /dev/null +++ b/examples/puredephasing.jl @@ -0,0 +1,126 @@ +#= + Example of a Pure Dephasing Model at finite or zero temperature with an hard cut-off Ohmic spectral density J(ω) = 2α(ω^s)/(ω_c^(s-1)) 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 + \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 + +import MPSDynamics: chaincoeffs_ohmic, puredephasingmpo + +const ∞ = Inf +#---------------------------- +# Physical parameters +#---------------------------- + +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(ω) + +#β = 100 # Thermalized environment + +β = ∞ # Case zero temperature T=0, β → ∞ + +#---------------------------- +# Ohmic spectral density +#---------------------------- + +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=true) # 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 + + +#----------------------- +# Simulation parameters +#----------------------- + +dt = 1.0 # time step + +tfinal = 10.0 # simulation time + +method = :TDVP1 # time-evolution method + +#method = :DTDVP # time-evolution method + +D = 2 # MPS bond dimension + +#--------------------------- +# MPO and initial state MPS +#--------------------------- + +H = puredephasingmpo(ω0, 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, prec=1E-4; + name = "pure dephasing model with temperature", + method = method, + obs = [ob1], + convobs = [ob1], + params = @LogParams(ω0, 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") + diff --git a/examples/sbm_Htimedependent.jl b/examples/sbm_Htimedependent.jl new file mode 100644 index 0000000..8d7896c --- /dev/null +++ b/examples/sbm_Htimedependent.jl @@ -0,0 +1,99 @@ +#= + Example of a time-dependent Hamitloninan on 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 + σ_x ϵ sin(ωdrive t) + 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 + +import MPSDynamics: disp + +#---------------------------- +# Physical parameters +#---------------------------- + +d = 4 # number of Fock states of the chain modes + +N = 60 # length of the chain + +α = 0.005 # coupling strength + +Δ = 0.0 # tunneling + +ω0 = 0.8 # TLS gap + +s = 1 # ohmicity + +cpars = chaincoeffs_ohmic(N, α, s) # chain parameters, i.e. on-site energies ϵ_i, hopping energies t_i, and system-chain coupling c_0 + +Trabi = 30.0 # Rabi period of the drive + +ϵ = 2*pi / Trabi # Intensity of the drive + +ωdrive = ω0 # Frequency of the drive + +Ndrive = 1 #Number of the site on which the drive is applied + +#----------------------- +# Simulation parameters +#----------------------- + +dt = 0.5 # time step + +tfinal = 100.0 # simulation time + +method = :TDVP1 # time-evolution method + +D = [6] # MPS bond dimension + +#--------------------------- +# MPO and initial state MPS +#--------------------------- + +timelist = collect(0:dt:tfinal) +numsteps = length(timelist)-1 + +Ht = [ϵ*sx*sin(ωdrive*tstep) for tstep in timelist] # Time-dependent Hamiltonian term + +H = spinbosonmpo(ω0, Δ, d, N, cpars) # MPO representation of the Hamiltonian + +ψ = unitcol(2,2) # Initial down-z system state + +A = productstatemps(physdims(H), state=[ψ, fill(unitcol(1,d), N)...]) # MPS representation of |ψ>|Vacuum> + +#--------------------------- +# Definition of observables +#--------------------------- + +ob1 = OneSiteObservable("sz", sz, 1) + +ob2 = OneSiteObservable("chain mode occupation", numb(d), (2,N+1)) + +ob3 = TwoSiteObservable("SXdisp", sx, disp(d), [1], collect(2:N+1)) + +#------------- +# Simulation +#------------ + +A, dat = runsim(dt, tfinal, A, H; + name = "Driving field on ohmic spin boson model", + method = method, + obs = [ob1], + convobs = [ob1], + params = @LogParams(N, d, α, Δ, ω0, s), + convparams = D, + timedep = true, # the Hamiltonian is time dependent + Ndrive = Ndrive, # the first site of the MPS/MPO (i.e. the system) is concerned + Htime = Ht, # list of time-dependent terms + verbose = false, + save = true, + plot = true, + ); + +#---------- +# Plots +#---------- + +plot(dat["data/times"], dat["data/sz"], label="Dmax = $(D...)", xlabel=L"t",ylabel=L"\sigma_z", title="") diff --git a/examples/sbm_zero_temperature.jl b/examples/sbm_zero_temperature.jl index bf8bdc9..6d0772f 100644 --- a/examples/sbm_zero_temperature.jl +++ b/examples/sbm_zero_temperature.jl @@ -1,10 +1,12 @@ -""" +#= 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 -""" + 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 + + Two variants of the one-site Time Dependent Variational Principal (TDVP) are presented for the time evolution of the quantum state. +=# using MPSDynamics, Plots, LaTeXStrings @@ -34,9 +36,13 @@ dt = 0.5 # time step tfinal = 30.0 # simulation time -method = :TDVP1 # time-evolution method +method = :TDVP1 # Regular one-site TDVP (fixed bond dimension) + +# method = :DTDVP # Adaptive one-site TDVP (dynamically updating bond dimension) + +convparams = [2,4,6] # MPS bond dimension (1TDVP) -D = [2,4,6] # MPS bond dimension +# convparams = [1e-2, 1e-3, 1e-4] # threshold value of the projection error (DTDVP) #--------------------------- # MPO and initial state MPS @@ -56,7 +62,7 @@ ob1 = OneSiteObservable("sz", sz, 1) ob2 = OneSiteObservable("chain mode occupation", numb(d), (2,N+1)) -ob3 = TwoSiteObservable("SXdisp", sx, disp(d), [1], collect(2:N+1)) +ob3 = TwoSiteObservable("SXdisp", sx, MPSDynamics.disp(d), [1], collect(2:N+1)) #------------- # Simulation @@ -68,8 +74,9 @@ A, dat = runsim(dt, tfinal, A, H; obs = [ob2,ob3], convobs = [ob1], params = @LogParams(N, d, α, Δ, ω0, s), - convparams = D, + convparams = convparams, verbose = false, + savebonddims = true, # this keyword argument enables the bond dimension at each time step to be saved when using DTDVP save = true, plot = true, ); @@ -78,6 +85,10 @@ A, dat = runsim(dt, tfinal, A, H; # Plots #---------- -plot(dat["data/times"], dat["convdata/sz"], label=["Dmax = 2" "Dmax = 4" "Dmax = 6"], xlabel=L"t",ylabel=L"\sigma_z") +method == :TDVP1 && plot(dat["data/times"], dat["convdata/sz"], label=["Dmax = 2" "Dmax = 4" "Dmax = 6"], xlabel=L"t",ylabel=L"\sigma_z") + +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") + +method == :DTDVP && heatmap(dat["data/times"], collect(0:N+1), dat["data/bonddims"], xlabel=L"t",ylabel="bond index") heatmap(dat["data/times"], collect(1:N), abs.(dat["data/SXdisp"][1,:,:]), xlabel=L"t",ylabel="chain mode") diff --git a/src/MPSDynamics.jl b/src/MPSDynamics.jl index e30c227..8fd1049 100644 --- a/src/MPSDynamics.jl +++ b/src/MPSDynamics.jl @@ -1,6 +1,6 @@ module MPSDynamics -using JLD, HDF5, Random, Dates, Plots, Printf, Distributed, LinearAlgebra, DelimitedFiles, KrylovKit, TensorOperations, GraphRecipes, SpecialFunctions, ITensors +using JLD, HDF5, Random, Dates, Plots, Printf, Distributed, LinearAlgebra, DelimitedFiles, KrylovKit, TensorOperations, GraphRecipes, SpecialFunctions, ITensors, Interpolations include("fundamentals.jl") include("reshape.jl") @@ -28,6 +28,7 @@ include("run_DTDVP.jl") include("run_A1TDVP.jl") include("chainA1TDVP.jl") include("switchmpo.jl") +include("finitetemperature.jl") """ runsim(dt, tmax, A, H; @@ -62,8 +63,6 @@ Propagate the MPS `A` with the MPO `H` up to time `tmax` in time steps of `dt`. * `name`: Used to describe the calculation. This name will appear in the log.txt file """ - - function runsim(dt, tmax, A, H; method=:TDVP1, machine=LocalMachine(), @@ -153,7 +152,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, tunnelingmpo export productstatemps, physdims, randmps, bonddims, elementmps @@ -173,5 +172,9 @@ export @LogParams export MPOtoVector, MPStoVector +export rhoreduced_2sites, rhoreduced_1site, protontransfermpo + +export chaincoeffs_finiteT, chaincoeffs_fermionic, fermionicspectraldensity_finiteT + end diff --git a/src/finitetemperature.jl b/src/finitetemperature.jl new file mode 100644 index 0000000..f4c4c5f --- /dev/null +++ b/src/finitetemperature.jl @@ -0,0 +1,212 @@ +using HDF5 +include("../ChainOhmT/quadohmT.jl") +include("../ChainOhmT/mcdis2.jl") + + +""" + chaincoeffs_finiteT(nummodes, β, ohmic=true; α, s, J, ωc=1, mc=4, mp=0, AB=nothing, iq=1, idelta=2, procedure=:Lanczos, Mmax=5000, save=true) + +Generate chain coefficients ``[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]`` for a harmonic bath at the inverse temperature β. + +By default a Ohmic spectral density ``J(ω) = 2αω_c (\\frac{ω}{ω_c})^s θ(ω-ω_c)`` is considered. +Users can provide their own spectral density. + +# Arguments +* nummodes: Number of bath modes +* β: inverse temperature +* ohmic: true if the spectral density is Ohmic, false if the user provides its own spectral density +* α: Kondo parameter of the Ohmic spectral density +* s: ohmicity +* J: user-provided spectral density. Should be a function f(x,i) where x is the frequency and i ∈ {1,...,mc} labels the intervals on which the SD is defined +* ωc: the maximum frequency of the Ohmic spectral density +* mc: the number of component intervals +* mp: the number of points in the discrete part of the measure (mp=0 if there is none) +* iq: a parameter to be set equal to 1, if the user provides his or her own quadrature routine, and different from 1 otherwise +* idelta: 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 +* procedure: choice between the Stieltjes and the Lanczos procedure +* AB: component intervals +* Mmax: maximum number of integration points +* save: if true the coefficients are saved +""" +function chaincoeffs_finiteT(nummodes, β, ohmic=true; α=1, s=1, J=nothing, ωc=1, mc=4, mp=0, AB=nothing, iq=1, idelta=2, procedure=:Lanczos, Mmax=5000, save=true) + + N = nummodes #Number of bath modes + + if AB==nothing + if mc==4 + AB = [[-Inf -ωc];[-ωc 0];[0 ωc];[ωc Inf]] + else + throw(ArgumentError("An interval AB with mc = $mc components should have been provided.")) + end + elseif length(AB) != mc + throw(ArgumentError("AB has a different number of intervals than mc = $mc.")) + end + + if ohmic==true + wf(x,i) = ohmicspectraldensity_finiteT(x,i,α,s,ωc,β) + elseif J==nothing + throw(ArgumentError("A spectral density should have been provided.")) + else + wf = J + end + + if procedure==:Lanczos # choice between the Stieltjes (irout = 1) and the Lanczos procedure (irout != 1) + irout = 2 + elseif procedure==:Stieltjes + irout = 1 + else + throw(ArgumentError("Procedure should be either Lanczos or Stieltjes.")) + end + + eps0=1e7*eps(Float64) + + jacerg = zeros(N,2) + + ab = 0. + ab, Mcap, kount, suc, uv = mcdis(N,eps0,quadfinT,Mmax,idelta,mc,AB,wf,mp,irout) + for m = 1:N-1 + jacerg[m,1] = ab[m,1] #site energy e + jacerg[m,2] = sqrt(ab[m+1,2]) #hopping parameter t + end + jacerg[N,1] = ab[N,1] + + eta = 0. + for i = 1:mc + xw = quadfinT(Mcap,i,uv,mc,AB,wf) + eta += sum(xw[:,2]) + end + jacerg[N,2] = sqrt(eta) # system-chain coupling c + + if save==true + # Write a HDF5 file + #curdir = @__DIR__ + dir = @__DIR__ + curdir = abspath(joinpath(dir, "../ChainOhmT")) + + if ohmic==true + Nstr=string(N) + astr=string(α) + sstr=string(s) + bstr=string(β) + # 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]) + + + else + Nstr = string(N) + wstr = string(ωc) + bstr = string(β) + # the "path" to the data inside of the h5 file is N -> ωc -> beta -> data (e, t or c) + + # Write onsite energies + h5write("$curdir/ohmicT/chaincoeffs.h5", string("/", Nstr, "/", wstr, "/", bstr, "/e"), jacerg[1:N,1]) + # Write hopping energies + h5write("$curdir/ohmicT/chaincoeffs.h5", string("/", Nstr, "/", wstr, "/", bstr,"/t"), jacerg[1:N-1,2]) + # Write coupling coefficient + h5write("$curdir/ohmicT/chaincoeffs.h5", string("/", Nstr, "/", wstr, "/", bstr,"/c"), jacerg[N,2]) + + end + end + + return [jacerg[:,1], jacerg[1:N-1,2],jacerg[N,2]] +end + +""" + chaincoeffs_fermionic(nummodes, β, chain; ϵ=x, ωc=1, mc=4, mp=0, AB=nothing, iq=1, idelta=2, procedure=:Lanczos, Mmax=5000, save=true) + +Generate chain coefficients ``[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]`` for a fermionic bath at the inverse temperature β. + +# Arguments +* nummodes: Number of bath modes +* β: inverse temperature +* chain: 1 if the chain modes are empty, 2 if the chain modes are filled +* ϵ: user-provided dispersion relation. Should be a function f(x) where x is the wavenumber +* J: user-provided spectral density. Should be a function f(x) where x is the wavenumber +* ωc: the maximum frequency allowwed in the spectral density +* mc: the number of component intervals +* mp: the number of points in the discrete part of the measure (mp=0 if there is none) +* iq: a parameter to be set equal to 1, if the user provides his or her own quadrature routine, and different from 1 otherwise +* idelta: 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 +* procedure: choice between the Stieltjes and the Lanczos procedure +* AB: component intervals +* Mmax: maximum number of integration points +* save: if true the coefficients are saved +""" + + +function chaincoeffs_fermionic(nummodes, β, chain; ϵ=nothing, J=nothing, ωc=1, mc=4, mp=0, AB=nothing, iq=1, idelta=2, procedure=:Lanczos, Mmax=5000, save=true) + + N = nummodes # Number of bath modes + + if AB==nothing + if mc==4 + AB = [[-Inf -ωc];[-ωc 0];[0 ωc];[ωc Inf]] + else + throw(ArgumentError("An interval AB with mc = $mc components should have been provided.")) + end + elseif length(AB) != mc + throw(ArgumentError("AB has a different number of intervals than mc = $mc.")) + end + + if ϵ==nothing + throw(ArgumentError("A dispersion relation should have been provided.")) + elseif J==nothing + throw(ArgumentError("The spectral density should be provided.")) + else + wf(x,i) = fermionicspectraldensity_finiteT(x, i , β, chain, ϵ, J) + end + + if procedure==:Lanczos # choice between the Stieltjes (irout = 1) and the Lanczos procedure (irout != 1) + irout = 2 + elseif procedure==:Stieltjes + irout = 1 + else + throw(ArgumentError("Procedure should be either Lanczos or Stieltjes.")) + end + + eps0=1e7*eps(Float64) + + jacerg = zeros(N,2) + + ab = 0. + ab, Mcap, kount, suc, uv = mcdis(N,eps0,quadfinT,Mmax,idelta,mc,AB,wf,mp,irout) + for m = 1:N-1 + jacerg[m,1] = ab[m,1] #site energy e + jacerg[m,2] = sqrt(ab[m+1,2]) #hopping parameter t + end + jacerg[N,1] = ab[N,1] + + eta = 0. + for i = 1:mc + xw = quadfinT(Mcap,i,uv,mc,AB,wf) + eta += sum(xw[:,2]) + end + jacerg[N,2] = sqrt(eta) # system-chain coupling c + + if save==true + # Write a HDF5 file + #curdir = @__DIR__ + dir = @__DIR__ + curdir = abspath(joinpath(dir, "../ChainOhmT")) + + Nstr=string(N) + cstr=string(chain) + bstr=string(β) + # the "path" to the data inside of the h5 file is beta -> alpha -> s -> data (e, t or c) + + # Write onsite energies + h5write("$curdir/fermionicT/chaincoeffs.h5", string("/", Nstr, "/", bstr, "/", cstr, "/e"), jacerg[1:N,1]) + # Write hopping energies + h5write("$curdir/fermionicT/chaincoeffs.h5", string("/", Nstr, "/", bstr, "/", cstr, "/t"), jacerg[1:N-1,2]) + # Write coupling coefficient + h5write("$curdir/fermionicT/chaincoeffs.h5", string("/", Nstr, "/", bstr, "/", cstr, "/c"), jacerg[N,2]) + end + + return [jacerg[:,1], jacerg[1:N-1,2],jacerg[N,2]] +end \ No newline at end of file diff --git a/src/fundamentals.jl b/src/fundamentals.jl index 87f686b..fcffbae 100644 --- a/src/fundamentals.jl +++ b/src/fundamentals.jl @@ -4,6 +4,7 @@ crea(d) = diagm(-1 => [sqrt(i) for i=1:d-1]) anih(d) = Matrix(crea(d)') numb(d) = crea(d)*anih(d) disp(d) = (1/sqrt(2))*(crea(d)+anih(d)) +disp(d,ωvib,m) = (1/(2*sqrt(m*ωvib/2)))*(crea(d)+anih(d)) mome(d) = (1/sqrt(2))*(im*(crea(d)-anih(d))) sx = [0. 1.; 1. 0.] sz = [1. 0.; 0. -1.] @@ -175,6 +176,125 @@ function measuremodes(adaga, e::Vector, t::Vector) return real.(diag(U' * adaga * U)) end +""" + measurecorrs(oper, , e::Vector, t::Vector) + + ### 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. + + ### 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``. + + ### 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. + +""" + +function measurecorrs(oper, e::Vector, t::Vector) + N = size(oper)[1] + hmat = diagm(0=>e[1:N], 1=>t[1:N-1], -1=>t[1:N-1]) + F = eigen(hmat) + U = F.vectors + return (U' * oper * U) +end + + +""" + cosineh(omega, bet) + + 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. + + # Returns + - `Float64`: The result of the modified cosine function. +""" + +function cosineh(omega, bet) + return 1/sqrt(1 - exp(-omega * (bet))) +end + +""" + sineh(omega, bet) + + 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. + + # Returns + - `Float64`: The result of the modified cosine function. +""" + +function sineh(omega, bet) + return 1/sqrt(-1 + exp(omega * float(bet))) +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. + + # 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. +""" + +function physical_occup(corr_constr, corr_destr, omega, occup, b, M) + x = range(-1, stop=1, length=M) + + # Ensure occup is a vector + occup_vector = vec(occup) + + occup_interp = LinearInterpolation(omega, occup_vector, extrapolation_bc=Line()) + corr_constr_interp = interpolate((omega, omega), abs.(corr_constr), Gridded(Linear())) + corr_destr_interp = interpolate((omega, omega), abs.(corr_destr), Gridded(Linear())) + + Mhalf = div(M, 2) + phys_occ = [] + + omega_min = minimum(omega) + omega_max = maximum(omega) + x_rescaled = (x .+ 1) .* (omega_max - omega_min) / 2 .+ omega_min + + for el in 1:Mhalf + ipos = Mhalf + el + ineg = Mhalf - el + 1 + occ = (cosineh(x_rescaled[ipos], b) * sineh(x_rescaled[ipos], b) * + (corr_destr_interp(x_rescaled[ineg], x_rescaled[ipos]) + corr_constr_interp(x_rescaled[ipos], x_rescaled[ineg])) + + cosineh(x_rescaled[ipos], b)^2 * occup_interp(x_rescaled[ipos]) + + sineh(x_rescaled[ipos], b)^2 * (1 + occup_interp(x_rescaled[ineg]))) + push!(phys_occ, occ) + end + + return phys_occ +end + + """ findchainlength(T, cparams...; eps=10^-6) diff --git a/src/measure.jl b/src/measure.jl index cc9b43b..92986c7 100644 --- a/src/measure.jl +++ b/src/measure.jl @@ -129,7 +129,7 @@ function measurempo(A::Vector, M::Vector, sites::Tuple{Int, Int}) for k=sites[1]:sites[2] F = updateleftenv(A[k], M[k], F) end - F = tensortrace(F, [1,2,1], [2]) + F = tensortrace([2], F, [1,2,1]) real(only(F)) end @@ -1027,3 +1027,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 ba42fea..bb28929 100644 --- a/src/models.jl +++ b/src/models.jl @@ -409,7 +409,7 @@ The spin is on site 1 of the MPS and the bath modes are to the right. This Hamiltonain is unitarily equivalent (before the truncation to `N` sites) to the spin-boson Hamiltonian defined by `` -H = \\frac{ω_0}{2}σ_z + Δσ_x + σ_x\\int_0^∞ dω\\sqrt{\\frac{J(ω)}{π}}(b_ω^\\dagger+b_ω) + \\int_0^∞ dω ωb_ω^\\dagger b_ω +H = \\frac{ω_0}{2}σ_z + Δσ_x + σ_x\\int_0^∞ dω\\sqrt{J(ω)}(b_ω^\\dagger+b_ω) + \\int_0^∞ dω ωb_ω^\\dagger b_ω ``. The chain parameters, supplied by `chainparams`=``[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]``, can be chosen to represent any arbitrary spectral density ``J(ω)`` at any temperature. @@ -486,8 +486,8 @@ end Generate chain coefficients ``[[ϵ_0,ϵ_1,...],[t_0,t_1,...],c_0]`` for an Harmonic bath at zero temperature with a power law spectral density given by: -soft cutoff: ``J(ω) = 2παω_c (\\frac{ω}{ω_c})^s \\exp(-ω/ω_c)`` \n -hard cutoff: ``J(ω) = 2παω_c (\\frac{ω}{ω_c})^s θ(ω-ω_c)`` +soft cutoff: ``J(ω) = 2αω_c (\\frac{ω}{ω_c})^s \\exp(-ω/ω_c)`` \n +hard cutoff: ``J(ω) = 2αω_c (\\frac{ω}{ω_c})^s θ(ω-ω_c)`` The coefficients parameterise the chain Hamiltonian @@ -498,8 +498,8 @@ H = H_S + c_0 A_S⊗B_0+\\sum_{i=0}^{N-1}t_i (b_{i+1}^\\dagger b_i +h.c.) + \\su which is unitarily equivalent (before the truncation to `N` sites) to `` -H = H_S + A_S⊗\\int_0^∞dω\\sqrt{\\frac{J(ω)}{π}}B_ω + \\int_0^∞dωωb_ω^\\dagger b_ω -`` +H = H_S + A_S⊗\\int_0^∞dω\\sqrt{J(ω)}B_ω + \\int_0^∞dωωb_ω^\\dagger b_ω +``. """ function chaincoeffs_ohmic(nummodes, α, s; ωc=1, soft=false) @@ -508,10 +508,16 @@ function chaincoeffs_ohmic(nummodes, α, s; ωc=1, soft=false) e = [ωc*(2n + 1 + s) for n in 0:(nummodes-1)] t = [ωc*sqrt((n + 1)*(n + s + 1)) for n in 0:(nummodes-2)] return [e, t, c0] - else - c0 = sqrt(2α/(s+1))*ωc - e = [(ωc/2)*(1 + (s^2)/((s+2n)*(2+s+2n))) for n in 0:(nummodes-1)] - t = [ωc*(1+n)*(1+s+n)/((s+2+2n)*(s+3+2n))*sqrt((3+s+2n)/(1+s+2n)) for n in 0:(nummodes-2)] + else + if s==0 + c0 = sqrt(2α)*ωc + e = fill(0,nummodes) + t = [ωc*(n+1)/(2n+1) for n in 0:(nummodes-2)] + else + c0 = sqrt(2α/(s+1))*ωc + e = [(ωc/2)*(1 + (s^2)/((s+2n)*(2+s+2n))) for n in 0:(nummodes-1)] + t = [ωc*(1+n)*(1+s+n)/((s+2+2n)*(s+3+2n))*sqrt((3+s+2n)/(1+s+2n)) for n in 0:(nummodes-2)] + end return [e, t, c0] end end @@ -699,7 +705,34 @@ function nearestneighbourmpo(tree_::Tree, h0, A, Ad = A') return TreeNetwork(tree, Ms) end +""" + puredephasingmpo(ΔE, dchain, Nchain, chainparams; tree=false) + + Generate MPO for a pure dephasing model, defined by the Hamiltonian + ``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 `` + + 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 +""" +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 """ tightbinding_mpo(N, ϵd, chainparams1, chainparams2) @@ -954,4 +987,409 @@ function interleaved_tightbinding_mpo(N, ϵd, chainparams1, chainparams2) M1fin[D,1,:,:] = u push!(W, M1fin) return W + +end + +""" + correlatedenvironmentmpo(R::Vector, Nm::Int, d::Int; chainparams, fnamecc::String, s=1, α=1, ωc=1, c_phonon=1, β="inf", issoft=false) + +Generate a MPO for a one-dimensional bosonic bath spatially correlated to a multi-component system + +`` +H_B + H_int = \\int_{-∞}^{+∞} dk ω_k b_k^\\dagger b_k + ∑_j \\int_{-∞}^{+∞}dk \\sqrt{J(k)}(A_j b_k e^{i k R_j} + h.c.) +``. + +The interactions between the system and the chain-mapped bath are long range, i.e. each site interacts with all the chain modes. The spectral density is assumed to be Ohmic ``J(ω) = 2αωc(ω/ωc)^s``. + +# Arguments + +* `R`: List of system's components positions +* `Nm`: Number of chain modes. The actual number of mode will be doubled to account for the left and right moving excitations. +* `d`: Local Hilbert space dimension of the bath modes +* `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. +* `fnamecc`: Path to a file containing pre-computed long-range coupling coefficient. If not provided, the coupling coefficients will be computed and stored. +* `s`: Ohmicity +* `α`: Kondo parameter +* `ωc`: Bath cut-off frequency +* `c_phonon`: Speed of sound in the bath +* `β`: Inverse temperature +* `issoft`: Is the cut-off of the Ohmic SD soft or hard? +""" +function correlatedenvironmentmpo(R::Vector, Nm::Int, d::Int; chainparams, fnamecc::String, s=1, α=1, ωc=1, c_phonon=1, β="inf", issoft=false) + + function polybeta(t::Float64, n::Int, a::Array, b::Array, temp::Array) + """ + polybeta recursively constructs the polynomials used to compute the coupling coefficients given the coefficients a and b + This function is useful when working at finite temperature (β != inf) + """ + if n==-1 + return 0 + elseif n==0 + if length(temp)>=2 + temp[2] = 1 + else + push!(temp,1) + end + + return 1 + elseif n==1 + pn = (t - a[n]) + if length(temp) == 2 + push!(temp,pn) + elseif length(temp) == 1 + push!(temp, 1) + push!(temp,pn) + end + + return pn + else + if length(temp) n+1 && temp[1] == t + pn = temp[n+2] + + return pn + else + temp = [t] + pn = (t - a[n])*polybeta(t,n-1,a,b,temp) - b[n-1]*polybeta(t,n-2,a,b,temp) #P_{n}(t) = (t-a_{n-1})P_{n-1} - b_{n-1}P_{n-2} + push!(temp, pn) + + return pn + end + end + end + + function SDOhmic(t) + """ + Bath Ohmic Spectral Density for zero temperature chain mapping of the bath + """ + if t==0 + return 0 + elseif t>-1 && t<1 + return 2*α*abs(t)*ωc + elseif abs(t)==1 + return 2*α*ωc + else + return 0 + end + end + + + function SDTOhmic(t) + """ + Bath Ohmic Spectral Density after the finite temperature chain mapping of the bath + """ + if t==0 + return 2*α/β + elseif t>-1 && t<1 + return α*t*ωc*(1+coth(β*t*ωc*0.5)) + elseif abs(t)==1 + return α*t*ωc*(1+coth(β*t*ωc*0.5)) + else + return 0 + end + end + + a_chain = chainparams[1] + b_chain = chainparams[2].^2 + + Norm = zeros(Nm) + + function γ(x::Int, n::Int, issoft::Bool; β="inf", temp=[1.]) + """ + Definition of the coupling coefficient between the site x and the mode n for a Ohmic spectral density with a hard cut-off (Jacobi Polynomials) or a soft cut-off Laguerre Polynomials + """ + if β=="inf" + if issoft==true + polynomial0(t) = sf_laguerre_n(n,s,t)*exp(-im*t*R[x]*ωc/c_phonon)*t^s*exp(-s) + return sqrt(2*α*gamma(n+s + 1)/gamma(n+1))*ωc*quadgk(polynomial0, 0, 1)[1] + else + polynomial(t) = jacobi(2*t-1,n-1, 0, s)*exp(-im*t*R[x]*ωc/c_phonon)*t^s + return sqrt(2*α*(2*(n-1) + s + 1))*ωc*quadgk(polynomial, 0, 1)[1] + end + elseif β!="inf" + polynomial1(t) = polybeta(t,n-1,a_chain,b_chain,[t]) + integrand(t) = polynomial1(t)*exp(im*t*R[x]*ωc/c_phonon)*SDTOhmic(t) + N2(t) = polynomial1(t)^2*SDTOhmic(t) + if Norm[n]==0 + Norm[n] = sqrt(quadgk(N2,-1,1)[1]) + end + return (ωc/Norm[n])*quadgk(integrand, -1, 1)[1] + + end + end + + # Construction of the MPO + W = Any[] # list of the MPO's tensors + + ### Construction of the bosonic bath MPO ### + e = chainparams[1] # list of the energy of each modes + t = chainparams[2] # list of the hopping energy between modes + u = unitmat(d) + # modes creation, anihilation and number operators + bd = crea(d) + b = anih(d) + n = numb(d) + + N = length(R) # Number of system's sites + + coupling_stored = zeros(ComplexF64,N,Nm) # just need NxNm because the two chains have the same coupling coeff up to complex conjugation + arestored = 1 # are the coupling coefficient stored + Nstored, Nmstored = N, Nm # number of stored coeff. + couplinglist = [] + try + couplinglist = readdlm(fnamecc,',',ComplexF64,'\n') + catch stored_error + if isa(stored_error, ArgumentError) + print("Coupling coefficient not found. They will be computed and stored.\n") + arestored = 0 + end + end + if couplinglist != [] + Nstored, Nmstored = size(couplinglist) + if Nmstored>=Nm && Nstored>=N + coupling_stored = couplinglist[1:N,1:Nm] + else + coupling_stored[1:min(N,Nstored),1:min(Nm,Nmstored)] = couplinglist[1:min(N,Nstored),1:min(Nm,Nmstored)] + print("Less coupling coefficient stored than needed. Available ones will be used and missing one will be computed and stored.\n") + end + end + + ## First chain MPO + D = 2*(N + 2) #Bond dimension + M = zeros(ComplexF64,D-2, D, d, d) + M[1,1,:,:] = M[D-2,D,:,:] = u + M[1,D,:,:] = e[1]*n + i = 2 + M[1,i,:,:] = t[1]*bd + M[1,i+1,:,:] = t[1]*b + + a = 0 #site counter + while i reshape(exp[i], size(exp[i])..., 1)) end + if reduceddensity + :Nrho in keys(kwargs) ? Nrho = kwargs[:Nrho] : Nrho = 1 + exprho = rhoreduced_1site(A0,Nrho) + push!(data, "Reduced ρ" => reshape(exprho, size(exprho)..., 1)) + end timed && (ttdvp = Vector{Float64}(undef, numsteps)) @@ -19,18 +25,27 @@ function run_1TDVP(dt, tmax, A, H, Dmax; obs=[], timed=false, kwargs...) for tstep=1:numsteps @printf("%i/%i, t = %.3f ", tstep, numsteps, times[tstep]) println() + if timedep + Ndrive = kwargs[:Ndrive] + Htime = kwargs[:Htime] + H0[Ndrive][1,1,:,:] = H[Ndrive][1,1,:,:] + Htime[tstep][:,:] + end if timed - val, t, bytes, gctime, memallocs = @timed tdvp1sweep!(dt, A0, H, F; kwargs...) + val, t, bytes, gctime, memallocs = @timed tdvp1sweep!(dt, A0, H0, F; kwargs...) println("\t","ΔT = ", t) A0, F = val ttdvp[tstep] = t else - A0, F = tdvp1sweep!(dt, A0, H, F; kwargs...) + A0, F = tdvp1sweep!(dt, A0, H0, F; kwargs...) end exp = measure(A0, obs; t=times[tstep]) 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,Nrho) + data["Reduced ρ"] = cat(data["Reduced ρ"], exprho; dims=ndims(exprho)+1) + end end timed && push!(data, "deltat"=>ttdvp) push!(data, "times" => times) diff --git a/src/run_2TDVP.jl b/src/run_2TDVP.jl index 2619b99..805ff87 100644 --- a/src/run_2TDVP.jl +++ b/src/run_2TDVP.jl @@ -1,5 +1,6 @@ -function run_2TDVP(dt, tmax, A, H, truncerr; obs=[], Dlim=50, savebonddims=false, timed=false, kwargs...) +function run_2TDVP(dt, tmax, A, H, truncerr; obs=[], Dlim=50, savebonddims=false, timed=false, reduceddensity=false, timedep=false, kwargs...) A0=deepcopy(A) + H0=deepcopy(H) data = Dict{String,Any}() numsteps = length(collect(0:dt:tmax))-1 @@ -11,6 +12,11 @@ function run_2TDVP(dt, tmax, A, H, truncerr; obs=[], Dlim=50, savebonddims=false for i=1:length(obs) push!(data, obs[i].name => reshape(exp[i], size(exp[i])..., 1)) end + if reduceddensity + :Nrho in keys(kwargs) ? Nrho = kwargs[:Nrho] : Nrho = 1 + exprho = rhoreduced_1site(A0,Nrho) + push!(data, "Reduced ρ" => reshape(exprho, size(exprho)..., 1)) + end bonds = bonddims(A0) savebonddims && push!(data, "bonddims" => reshape([bonds...], length(bonds), 1)) @@ -22,19 +28,28 @@ function run_2TDVP(dt, tmax, A, H, truncerr; obs=[], Dlim=50, savebonddims=false maxbond = max(bonds...) @printf("%i/%i, t = %.3f, Dmax = %i ", tstep, numsteps, times[tstep], maxbond) println() + if timedep + Ndrive = kwargs[:Ndrive] + Htime = kwargs[:Htime] + H0[Ndrive][1,1,:,:] = H[Ndrive][1,1,:,:] + Htime[tstep][:,:] + end if timed - val, t, bytes, gctime, memallocs = @timed tdvp2sweep!(dt, A0, H, F; truncerr=truncerr, truncdim=Dlim, kwargs...) + val, t, bytes, gctime, memallocs = @timed tdvp2sweep!(dt, A0, H0, F; truncerr=truncerr, truncdim=Dlim, kwargs...) println("\t","ΔT = ", t) A0, F = val ttdvp[tstep] = t else - A0, F = tdvp2sweep!(dt, A0, H, F; truncerr=truncerr, truncdim=Dlim, kwargs...) + A0, F = tdvp2sweep!(dt, A0, H0, F; truncerr=truncerr, truncdim=Dlim, kwargs...) end bonds = bonddims(A0) exp = measure(A0, obs; t=times[tstep]) 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,Nrho) + data["Reduced ρ"] = cat(data["Reduced ρ"], exprho; dims=ndims(exprho)+1) + end if savebonddims data["bonddims"] = cat(data["bonddims"], [bonds...], dims=2) end diff --git a/src/run_DTDVP.jl b/src/run_DTDVP.jl index c39d34e..928abe4 100644 --- a/src/run_DTDVP.jl +++ b/src/run_DTDVP.jl @@ -1,5 +1,6 @@ -function run_DTDVP(dt, tmax, A, H, prec; obs=[], effects=false, error=false, timed=false, savebonddims=false, Dplusmax=nothing, Dlim=50, kwargs...) +function run_DTDVP(dt, tmax, A, H, prec; obs=[], effects=false, error=false, timed=false, savebonddims=false, Dplusmax=nothing, Dlim=50, reduceddensity=false, timedep=false, kwargs...) A0=deepcopy(A) + H0=deepcopy(H) data = Dict{String,Any}() numsteps = length(collect(0:dt:tmax))-1 @@ -11,6 +12,11 @@ function run_DTDVP(dt, tmax, A, H, prec; obs=[], effects=false, error=false, tim for i=1:length(obs) push!(data, obs[i].name => reshape(exp[i], size(exp[i])..., 1)) end + if reduceddensity + :Nrho in keys(kwargs) ? Nrho = kwargs[:Nrho] : Nrho = 1 + exprho = rhoreduced_1site(A0,Nrho) + push!(data, "Reduced ρ" => reshape(exprho, size(exprho)..., 1)) + end bonds = bonddims(A0) savebonddims && push!(data, "bonddims" => reshape([bonds...], length(bonds), 1)) @@ -25,7 +31,12 @@ function run_DTDVP(dt, tmax, A, H, prec; obs=[], effects=false, error=false, tim for tstep=1:numsteps maxbond = max(bonds...) @printf("%i/%i, t = %.3f, Dmax = %i \n", tstep, numsteps, times[tstep], maxbond) - A0, Afull, F, info = tdvp1sweep_dynamic!(dt, A0, H, Afull, F; + if timedep + Ndrive = kwargs[:Ndrive] + Htime = kwargs[:Htime] + H0[Ndrive][1,1,:,:] = H[Ndrive][1,1,:,:] + Htime[tstep][:,:] + end + A0, Afull, F, info = tdvp1sweep_dynamic!(dt, A0, H0, Afull, F; obs=obs, prec=prec, Dlim=Dlim, @@ -46,6 +57,10 @@ function run_DTDVP(dt, tmax, A, H, prec; obs=[], effects=false, error=false, tim 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,Nrho) + data["Reduced ρ"] = cat(data["Reduced ρ"], exprho; dims=ndims(exprho)+1) + end end if savebonddims data["bonddims"] = cat(data["bonddims"], bonds, dims=2) @@ -55,6 +70,11 @@ function run_DTDVP(dt, tmax, A, H, prec; obs=[], effects=false, error=false, tim 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,Nrho) + data["Reduced ρ"] = cat(data["Reduced ρ"], exprho; dims=ndims(exprho)+1) + end + if effects efftarray = efft[1] diff --git a/src/treeTDVP.jl b/src/treeTDVP.jl index 26c87dd..cddc2e6 100644 --- a/src/treeTDVP.jl +++ b/src/treeTDVP.jl @@ -50,8 +50,9 @@ function mpsrightnorm!(net::TreeNetwork, id::Int) IC=collect(1:nc+2) IA=collect(1:nc+2) IC[i+1]=-1 + #net[id] = tensorcontract(net[id], IA, C, [i+1,-1], IC) #Old tensoroperation version + net[id] = tensorcontract(IC,net[id], IA, C, [i+1,-1]) - net[id] = tensorcontract(net[id], IA, C, [i+1,-1], IC) end end @@ -405,7 +406,8 @@ function tdvp1sweep_lc!(dt, A::TreeNetwork, M::TreeNetwork, lc::TreeLightCone, F IA = collect(1:ngc+2) IB = collect(1:ngc+2) IA[1] = -1 - A[child] = tensorcontract(A[child], IA, C, [1,-1], IB) + #A[child] = tensorcontract(A[child], IA, C, [1,-1], IB) #old tensoroperations + A[child] = tensorcontract(IB, A[child], IA, C, [1,-1]) #(OC is now on child) #evolve child forward one full time step @@ -427,7 +429,8 @@ function tdvp1sweep_lc!(dt, A::TreeNetwork, M::TreeNetwork, lc::TreeLightCone, F IA = collect(1:nc+2) IB = collect(1:nc+2) IA[i+1] = -1 - AC = tensorcontract(AL, IA, C, [i+1,-1], IB) + #AC = tensorcontract(AL, IA, C, [i+1,-1], IB) #old tensoroperations + AC = tensorcontract(IB, AL, IA, C, [i+1,-1]) #(OC is now on headnode) end @@ -496,7 +499,8 @@ function tdvp1sweep_lc!(dt, A::TreeNetwork, M::TreeNetwork, lc::TreeLightCone, F IA = collect(1:ngc+2) IB = collect(1:ngc+2) IA[1] = -1 - A[child] = tensorcontract(A[child], IA, C, [1,-1], IB) + #A[child] = tensorcontract(A[child], IA, C, [1,-1], IB) + A[child] = tensorcontract(IB, A[child], IA, C, [1,-1]) #(OC is now on child) #evolve child forward one full time step @@ -518,7 +522,8 @@ function tdvp1sweep_lc!(dt, A::TreeNetwork, M::TreeNetwork, lc::TreeLightCone, F IA = collect(1:nc+2) IB = collect(1:nc+2) IA[i+1] = -1 - AC = tensorcontract(AL, IA, C, [i+1,-1], IB) + #AC = tensorcontract(AL, IA, C, [i+1,-1], IB) + AC = tensorcontract(IB, AL, IA, C, [i+1,-1]) #(OC is now on node) end