Merge branch 'master' into doc-writing

ChainOhmT/chainOhmT.jl

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

ChainOhmT/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
 

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

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

ChainOhmT/quadohmT.jl

@@ -1,73 +1,96 @@
-function quadfinT(N,i,uv)
+function quadfinT(N,i,uv,mc,AB,wf)
 
     if (i>1 && i<mc)
-        xw = tr(uv,i)
+        xw = trr(uv,i,AB,wf)
         return xw
     end
 
     if mc==1
       if (AB[i,1]!=-Inf)&&(AB[i,2]!=Inf)
-        xw = tr(uv,i)
+        xw = trr(uv,i,AB,wf)
         return xw
       elseif AB[i,1]!=-Inf
-          xw = rtr(uv,i)
+          xw = rtr(uv,i,AB,wf)
           return xw
       elseif AB[i,2]!=Inf
-          xw = ltr(uv,i)
+          xw = ltr(uv,i,AB,wf)
           return xw
       else
-          xw = str(uv,i)
+          xw = str(uv,i,wf)
           return xw
       end
     else
       if ((i==1 && AB[i,1]!=-Inf)||(i==mc && AB[i,2]!=Inf))
-        xw = tr(uv,i)
+        xw = trr(uv,i,AB,wf)
         return xw
       elseif i==1
-        xw = ltr(uv,i)
+        xw = ltr(uv,i,AB,wf)
         return xw
       end
     end
 
-    xw = rtr(uv,mc)
+    xw = rtr(uv,mc,AB,wf)
     return xw
 end
 
-function tr(t,i)
+function trr(t,i,AB,wf)
     s = Array{Float64}(undef, size(t))
     s[:,1] = ((AB[i,2]-AB[i,1])*t[:,1] .+ AB[i,2] .+ AB[i,1])./2
     s[:,2] = (AB[i,2]-AB[i,1]).*t[:,2].*wf(s[:,1],i)./2
     return s
 end
 
-function str(t,i)
+function str(t,i,wf)
     s = Array{Float64}(undef, size(t))
     s[:,1] = t[:,1]./(1-t[:,1].^2)
     s[:,2] = t[:,2].*wf(s[:,1],i).*(1+t[:,1].^2)./((1-t[:,1].^2).^2)
     return s
 end
 
-function rtr(t,i)
+function rtr(t,i,AB,wf)
     s = Array{Float64}(undef, size(t))
     s[:,1] = AB[i,1] .+ (t[:,1] .+ 1)./(-t[:,1] .+ 1)
     s[:,2] = 2*t[:,2].*wf(s[:,1],i)./((1 .-t[:,1]).^2)
     return s
 end
 
-function ltr(t,i)
+function ltr(t,i,AB,wf)
     s = Array{Float64}(undef, size(t))
     s[:,1] = -(-t[:,1].+1)./(t[:,1].+1) .+ AB[i,2]
     s[:,2] = 2*t[:,2].*wf(s[:,1],i)./((t[:,1] .+ 1).^2)
     return s
 end
 
-function wf(x,i)
+function ohmicspectraldensity_finiteT(x,i,α,s,ωc,β)
     if i==1
         y = 0
     elseif i==2
-        y = -pi*a*abs.(x) .* (coth.((beta/2).*x) .+ 1) #.* exp(-abs(x)/wc)
+        y = -2 .*( α*abs.(x).^s ./ ωc^(s-1)) .* (coth.((β/2).*x) .+ 1)./2 #.* exp(-abs(x)/wc)
     elseif i==3
-        y = pi*a*abs.(x) .* (coth.((beta/2).*x) .+ 1) #.* exp(-abs(x)/wc)
+        y = 2 .*( α*abs.(x).^s ./ ωc^(s-1)) .* (coth.((β/2).*x) .+ 1)./2 #.* exp(-abs(x)/wc)
+    elseif i==4
+        y = 0
+    end
+    return y
+end
+
+"""
+    fermionicspectraldensity_finiteT(x, i, β, chain; ϵ=x)
+
+    Fermionic spectral density at finite temperature β. A function f(x,i) where x is the frequency 
+    and i ∈ {1,...,mc} labels the intervals on which the SD is defined.
+
+"""
+
+function fermionicspectraldensity_finiteT(x, i, β, chain, ϵ, J)
+    if i==1
+        y = 0
+    elseif i==2 || i==3
+        if chain==1
+            y = (J.(x) .* sqrt.(1. ./(exp.(-β .* ϵ.(x)) .+ 1))).^2
+        elseif chain==2
+            y = (J.(x) .* sqrt.(1. ./(exp.(β .* ϵ.(x)) .+ 1))).^2
+        end
     elseif i==4
         y = 0
     end

ChainOhmT/stieltjes.jl

@@ -9,8 +9,8 @@
     coefficients in the second column, of the nx2 array ab.
 """
 function stieltjes(N,xw) #return ab
-    tiny = 10*realmin
-    huge = .1*realmax
+    #tiny = 10*realmin
+    #huge = .1*realmax
 
     ## Remove data with zero weights ##
 
@@ -19,15 +19,16 @@ function stieltjes(N,xw) #return ab
     xw = xw[I,:]
     index = findfirst(x->x!=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)
 
     if (N<=0||N>Ncap)
         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)<tiny
-        error("impending underflow in stieltjes for k=$(k)")
-      end
+      p2 = (xw[:,1] .- ab[k,1]).*p1 .- ab[k,2].*p0
+      s1 = (transpose(xw[:,2])*(p2.^2))[1]
+      s2 = (transpose(xw[:,1])*(xw[:,2].*(p2.^2)))[1]
+      # if (max(abs(p2))>huge)||(abs(s2)>huge)
+      #   error("impending overflow in stieltjes for k=$(k)")
+      # end
+      # if abs(s1)<tiny
+      #   error("impending underflow in stieltjes for k=$(k)")
+      # end
       ab[k+1,1] = s2/s1
       ab[k+1,2] = s1/s0
       s0 = s1