jobs/CMPO_TFIsing.jl

 using LinearAlgebra; BLAS.set_num_threads(Threads.nthreads())
 using TimerOutputs, Dates, ArgParse
 using DelimitedFiles, HDF5, Printf
-using JuliaCMPO
+using JuliaCMPO, FiniteTLanczos
 
 println("JULIA_NUM_THREADS = ", Threads.nthreads())
 println("OPENBLAS_NUM_THREADS = ", BLAS.get_num_threads())
 
 const to = TimerOutput()
 const Start_Time = Dates.format(now(), "yyyy-mm-dd HH:MM:SS")
-settings = ArgParseSettings(prog="CMPO: TFIsing model"
+settings = ArgParseSettings(prog="CMPO: 1D TFIsing model"
 )
 @add_arg_table! settings begin
     "--J"
@@ -47,6 +47,14 @@ settings = ArgParseSettings(prog="CMPO: TFIsing model"
         arg_type = String
         default = "/data/sliang/JuliaCMPO/TFIsing"
         help = "result folder"
+    "--tag"
+        arg_type = String
+        default = Dates.format(now(), "yyyy-mm-dd")
+        help = "date tag"
+    "--processor"
+        arg_type = Int64
+        default = 0
+        help = "processor used: 0 for CPU and 1 for GPU"
 end
 
 parsed_args = parse_args(settings; as_symbols=true)
@@ -61,6 +69,9 @@ const β0 = parsed_args[:init]
 const bondD = parsed_args[:bondD]
 const Continue = parsed_args[:Continue]
 const ResultFolder = parsed_args[:ResultFolder]
+const tag = parsed_args[:tag]
+const processor = Processor(parsed_args[:processor])
+
 
 const wid = 1
 
@@ -71,14 +82,20 @@ isdir(ModelResultFolder) || mkdir(ModelResultFolder)
 #CMPO
 model = TFIsing(J, Γ)
 
+#trace_estimator = nothing
+#trace_estimator = TraceEstimator(Full_ED, FTLMOptions(processor = processor))
+trace_estimator = TraceEstimator(FullSampling_FTLM, FTLMOptions(Nk = 100, processor = processor))
+
+trace_estimator === nothing ? estimator_name = "nothing" : estimator_name = string(trace_estimator.estimator)
+
 if β0 == 0
     ψ0 = nothing
 else
-    init_path = @sprintf "%s/bondD_%02i_CMPS/beta_%.2f.hdf5" ModelResultFolder bondD β0
+    init_path = @sprintf "%s/bondD_%02i_CMPS_%s_%s/beta_%.2f.hdf5" ModelResultFolder bondD estimator_name tag β0
     ψ0 = readCMPS(init_path)
 end
 
-EngFile = @sprintf "%s/Obsv_FECvS_bondD_%02i.txt" ModelResultFolder bondD
+EngFile = @sprintf "%s/Obsv_FECvS_bondD_%02i_%s_%s.txt" ModelResultFolder bondD estimator_name tag
 if Continue == false
     open(EngFile,"w") do file  
         write(file, "  β          free_energy           energy              specific_heat            entropy      \n")
@@ -87,10 +104,16 @@ if Continue == false
 end
 
 βlist = [i for i in range(bi, bf, step=bstep)]
+
 for b = 1:length(βlist)
     β = βlist[b]
     @timeit to "evaluate" begin
-        res = evaluate(model, bondD, β, ModelResultFolder, init = ψ0)
+        evaluate_options = EvaluateOptions(trace_estimator=trace_estimator,
+                            init = ψ0,
+                            tag = tag,
+                            processor = processor)
+        res = JuliaCMPO.evaluate(model, bondD, β, ModelResultFolder, 
+                       options = evaluate_options)
     end
 
     open(EngFile,"a") do file  

jobs/CMPO_TFIsing_2D_helical.jl

@@ -43,6 +43,14 @@ settings = ArgParseSettings(prog="CMPO code for TFIsing_2D_helical"
         arg_type = String
         default = "/data/sliang/JuliaCMPO/TFIsing_2D_helical_unexpanded"
         help = "result folder"
+    "--tag"
+        arg_type = String
+        default = Dates.format(now(), "yyyy-mm-dd")
+        help = "date tag"
+    "--processor"
+        arg_type = Processor
+        default = CPU
+        help = "processor used"
 end
 parsed_args = parse_args(settings; as_symbols=true)
 print(parsed_args,"\n")
@@ -53,7 +61,12 @@ const width = parsed_args[:width]
 const β = parsed_args[:beta]
 const bondD = parsed_args[:bondD]
 const max_pow_step = parsed_args[:max_pow_step]
+const tag = parsed_args[:tag]
+const processor = parsed_args[:processor]
+
 Continue = parsed_args[:Continue]
+if Continue ≥ max_pow_step Continue = true end
+
 #Creat ResultFolders if there is none
 const ResultFolder = parsed_args[:ResultFolder]
 isdir(ResultFolder) || mkdir(ResultFolder)
@@ -63,11 +76,17 @@ isdir(ModelResultFolder) || mkdir(ModelResultFolder)
 
 #CMPO
 model = TFIsing_2D_helical(J, Γ, width, expand=false)
-if Continue ≥ max_pow_step Continue = true end
+
 @timeit to "evaluate" begin
-    res = evaluate(model, bondD, β, ModelResultFolder, 
-                        hermitian = false, 
-                        max_pow_step = max_pow_step, Continue = Continue)
+    evaluate_options = EvaluateOptions(EvaluateOptions(),
+                            init = ψ0,
+                            hermitian = false,
+                            max_pow_step = max_pow_step,
+                            Continue = Continue,
+                            tag = tag,
+                            processor = processor)
+    res = JuliaCMPO.evaluate(model, bondD, β, ModelResultFolder, 
+                            options = evaluate_options)
 end
 
 const End_Time = Dates.format(now(), "yyyy-mm-dd HH:MM:SS")

jobs/CMPO_XXZ.jl

@@ -55,10 +55,10 @@ settings = ArgParseSettings(prog="CMPO code for XXZ model"
         arg_type = String
         default = Dates.format(now(), "yyyy-mm-dd")
         help = "date tag"
-    "--device"
-        arg_type = Symbol
-        default = :cpu
-        help = "device used"
+    "--processor"
+        arg_type = Processor
+        default = CPU
+        help = "processor used"
 end
 parsed_args = parse_args(settings; as_symbols=true)
 print(parsed_args,"\n")
@@ -73,7 +73,7 @@ const bondD = parsed_args[:bondD]
 const init = parsed_args[:init]
 const Continue = parsed_args[:Continue]
 const tag = parsed_args[:tag]
-const device = parsed_args[:device]
+const processor = parsed_args[:processor]
 
 const wid = 1
 const ResultFolder = parsed_args[:ResultFolder]
@@ -84,7 +84,6 @@ isdir(ModelResultFolder) || mkdir(ModelResultFolder)
 
 #CMPO
 model = XXZmodel(Jz/Jxy)
-#group = 2  #AFM XXZ
 
 if β0 == 0
     ψ0 = nothing
@@ -102,15 +101,17 @@ if Continue == false
 end
 
 βlist = [i for i in range(bi, bf, step=bstep)]
-device == :cpu ? solver = cpu_solver : solver = gpu_solver
+
 for b = 1:length(βlist)
     β = βlist[b]
     @timeit to "evaluate" begin
-        res = evaluate(model, bondD, β, ModelResultFolder, 
-                            init = ψ0, 
-                            group = group, 
-                            tag=tag,
-                            solver = solver)
+        evaluate_options = EvaluateOptions(EvaluateOptions(),
+                            init = ψ0,
+                            group = group,
+                            tag = tag,
+                            processor = processor)
+        res = JuliaCMPO.evaluate(model, bondD, β, ModelResultFolder, 
+                       options = evaluate_options)
     end
 
     open(EngFile,"a") do file  

jobs/CMPO_XXZ_2D_helical.jl

-using LinearAlgebra; BLAS.set_num_threads(Threads.nthreads())
+using LinearAlgebra; BLAS.set_num_threads(8)
 using TimerOutputs, Dates, ArgParse
 using DelimitedFiles, HDF5, Printf
-using JuliaCMPO
+using JuliaCMPO, FiniteTLanczos
 
 println("JULIA_NUM_THREADS = ", Threads.nthreads())
 println("OPENBLAS_NUM_THREADS = ", BLAS.get_num_threads())
 
 const to = TimerOutput()
 const Start_Time = Dates.format(now(), "yyyy-mm-dd HH:MM:SS")
-settings = ArgParseSettings(prog="CMPO code for XXZ model"
-)
+settings = ArgParseSettings(prog="CMPO code for XXZ model")
+
 @add_arg_table! settings begin
     "--Jxy"
         arg_type = Float64
@@ -55,10 +55,10 @@ settings = ArgParseSettings(prog="CMPO code for XXZ model"
         arg_type = String
         default = Dates.format(now(), "yyyy-mm-dd")
         help = "date tag"
-    "--device"
-        arg_type = Symbol
-        default = :CPU
-        help = "device used"
+    "--processor"
+        arg_type = Int64
+        default = 0
+        help = "processor used: 0 for CPU and 1 for GPU"
 end
 parsed_args = parse_args(settings; as_symbols=true)
 print(parsed_args,"\n")
@@ -72,9 +72,11 @@ const β = parsed_args[:beta]
 const bondD = parsed_args[:bondD]
 const max_pow_step = parsed_args[:max_pow_step]
 const tag = parsed_args[:tag]
-const device = parsed_args[:device]
+const processor = Processor(parsed_args[:processor])
 
 Continue = parsed_args[:Continue]
+if Continue ≥ max_pow_step Continue = true end
+
 #Creat ResultFolders if there is none
 const ResultFolder = parsed_args[:ResultFolder]
 isdir(ResultFolder) || mkdir(ResultFolder)
@@ -85,18 +87,25 @@ isdir(ModelResultFolder) || mkdir(ModelResultFolder)
 #CMPO
 model = XXZmodel_2D_helical(Jz/Jxy, width, expand = expand)
 
-if Continue ≥ max_pow_step Continue = true end
-device == :CPU ? solver = cpu_solver : solver = gpu_solver
+trace_estimator = nothing
+#trace_estimator = TraceEstimator(simple_FTLM, FTLMOptions(processor = processor))
+#trace_estimator = TraceEstimator(orthogonalized_FTLM, FTLMOptions(processor = processor))
+
+trace_estimator === nothing ? estimator_name = "nothing" : estimator_name = string(trace_estimator.estimator)
+
 @timeit to "evaluate" begin
-    res = evaluate(model, bondD, β, ModelResultFolder, 
-                        hermitian = false, 
-                        max_pow_step = max_pow_step, 
-                        group = group,
-                        Continue = Continue,
-                        tag = tag,
-                        solver = solver)
+    evaluate_options = EvaluateOptions(trace_estimator = trace_estimator,
+                            hermitian = false,
+                            group = group,
+                            max_pow_step = max_pow_step,
+                            Continue = Continue,
+                            tag = tag,
+                            processor = processor)
+    res = JuliaCMPO.evaluate(model, bondD, β, ModelResultFolder, 
+                            options = evaluate_options)
 end
 
+
 const End_Time = Dates.format(now(), "yyyy-mm-dd HH:MM:SS")
 const Running_TimeTable = string(to)
 @show Start_Time

jobs/bright90.jl

@@ -2,14 +2,14 @@ using Printf, DelimitedFiles
 using Random; Random.seed!()
 
 @enum Partitian a100 p100 v100 titanv
-@enum Device CPU GPU
+@enum Processor CPU GPU
 
 """
     submitJob: Prepare a jobfile
 """
 function submitJob(env, prog, args, jobname;
                   partitian::Partitian = a100,
-                  device::Device = CPU,
+                  processor::Processor = CPU,
                   gpu_memory::Int = 0,
                   cpu_per_task::Integer = 4,
                   Run=false, 
@@ -30,7 +30,7 @@ function submitJob(env, prog, args, jobname;
     #SBATCH --error=%s 
     """ partitian runtime cpu_per_task jobname log_file_path log_file_path
     
-    if device == GPU
+    if processor == GPU
         if partitian == a100 && gpu_memory != 0
             job *= """
             #SBATCH --gres=gpu:A100_$(gpu_memory)G:1
@@ -89,7 +89,7 @@ function submitJob(env, prog, args, jobname;
 
     if Run 
         run(`sbatch $(jobfile)`) 
-        sleep(2.0)
+        sleep(1.0)
         if return_id
             jobid = readdlm(jobname*"_id.txt", Int64)[1]
             rm(jobname*"_id.txt")

jobs/jobsCMPO_TFIsing.jl

-using Printf
+using Printf, Dates
 include("/home/sliang/JuliaCode/JuliaCMPO/jobs/bright90.jl")
 
 phys_model = "TFIsing"
 env = "/home/sliang/JuliaCode/JuliaCMPO"
 prog = env * "/jobs/CMPO_$(phys_model).jl"
 
-bi = 11.3
-bf = 40.0
+processor = CPU
+machine = v100
+gpu_memory = 80
+
+logtag = Dates.format(now(), "yyyy-mm-dd")*"-"*"$(processor)"
+Wait = nothing
+cpu_per_task = 8
+
+tag = Dates.format(now(), "yyyy-mm-dd")*"-"*"$(processor)"
+
+bi = 1.0
+bf = 5.0
 bstep = 0.1
-init = bi - bstep
-Continue = true
+init = 0.0
+Continue = false
 
 Jlist = [1.0]
 Γlist = [1.0]
-init = 11.2
-bondDlist = [16]
-
+bondDlist = [8]
 
 
 #CREAT LOG FOLDER
-logdir = "/data/sliang/log/JuliaCMPO"
+logdir = "/data/sliang/JuliaCMPO/log"
 isdir(logdir) || mkdir(logdir)
 
 for J in Jlist, Γ in Γlist, bondD in bondDlist
@@ -30,15 +38,24 @@ for J in Jlist, Γ in Γlist, bondD in bondDlist
                     "bf"=>bf,
                     "bstep"=>bstep,
                     "init"=>init,
-                    "Continue"=>Continue)
+                    "Continue"=>Continue,
+                    "tag"=>tag,
+                    "processor"=> Int(processor)
+                    )
     jobname = logdir * "/" * phys_model
     isdir(jobname) || mkdir(jobname)
     jobname = @sprintf "%s/J_%.2f_G_%.2f_wid_01" jobname J Γ
     isdir(jobname) || mkdir(jobname)
-    jobname = @sprintf "%s/bondD_%02i" jobname bondD
+    jobname = @sprintf "%s/bondD_%02i_%s" jobname bondD logtag
     isdir(jobname) || mkdir(jobname)
-    jobname = @sprintf "%s/bi_%.2f_bf_%.2f_bstep_%.2f" jobname bi bf bstep
+    timetag = Dates.format(now(), "HH-MM-SS")
+    jobname = @sprintf "%s/bi_%.2f_bf_%.2f_bstep_%.2f_%s" jobname bi bf bstep timetag
 
     jobid = submitJob(env, prog, args, jobname, 
-                            Run = true)
+                            partitian = machine,
+                            processor = processor,
+                            gpu_memory = gpu_memory,
+                            cpu_per_task = cpu_per_task,
+                            Run = true, 
+                            Wait = Wait)
 end

jobs/jobsCMPO_TFIsing_2D_helical.jl

@@ -5,7 +5,16 @@ phys_model = "TFIsing_2D_helical"
 env = "/home/sliang/JuliaCode/JuliaCMPO"
 prog = env * "/jobs/CMPO_$(phys_model).jl"
 
+processor = GPU
+machine = a100
+gpu_memory = 80
+logtag = Dates.format(now(), "yyyy-mm-dd")*"-"*"$(processor)"
+
 Wait = nothing
+cpu_per_task = 8
+
+#tag = "2022-06-21"*"-"*processor
+tag = Dates.format(now(), "yyyy-mm-dd")*"-"*"$(processor)"
 
 βlist = [0.1]
 Jlist = [1.0]
@@ -27,7 +36,9 @@ for bondD in bondDlist
                     "width"=>width,
                     "beta"=>β,
                     "max_pow_step"=>max_pow_step,
-                    "Continue"=>Continue
+                    "Continue"=>Continue,
+                    "tag"=>tag,
+                    "processor"=>processor
                     )
         jobname = logdir * "/" * phys_model
         isdir(jobname) || mkdir(jobname)
@@ -37,6 +48,12 @@ for bondD in bondDlist
         isdir(jobname) || mkdir(jobname)
         jobname = @sprintf "%s/beta_%.2f" jobname β
 
-        jobid = submitJob(env, prog, args, jobname, Run = true, Wait = Wait)
+        jobid = submitJob(env, prog, args, jobname, 
+                            partitian = machine,
+                            processor = processor,
+                            gpu_memory = gpu_memory,
+                            cpu_per_task = cpu_per_task,
+                            Run = true, 
+                            Wait = Wait)
     end   
 end      

jobs/jobsCMPO_XXZ.jl

@@ -5,6 +5,17 @@ phys_model = "XXZ"
 env = "/home/sliang/JuliaCode/JuliaCMPO"
 prog = env * "/jobs/CMPO_$(phys_model).jl"
 
+processor = GPU
+machine = a100
+gpu_memory = 80
+logtag = Dates.format(now(), "yyyy-mm-dd")*"-"*"$(processor)"
+
+Wait = nothing
+cpu_per_task = 8
+
+#tag = "2022-06-21"*"-"*processor
+tag = Dates.format(now(), "yyyy-mm-dd")*"-"*"$(processor)"
+
 bi = 8.5
 bf = 20.0
 bstep = 0.1
@@ -29,7 +40,10 @@ for Jz in Jzlist, Jxy in Jxylist, bondD in bondDlist
                 "bf"=>bf,
                 "bstep"=>bstep,
                 "init"=>init,
-                "Continue"=>Continue)
+                "Continue"=>Continue,
+                "tag"=>tag,
+                "processor"=> processor
+                )
     jobname = logdir * "/" * phys_model
     isdir(jobname) || mkdir(jobname)
     jobname = @sprintf "%s/Jz_%.2f_Jxy_%.2f_wid_01" jobname Jz Jxy
@@ -38,5 +52,11 @@ for Jz in Jzlist, Jxy in Jxylist, bondD in bondDlist
     isdir(jobname) || mkdir(jobname)
     jobname = @sprintf "%s/bi_%.2f_bf_%.2f_bstep_%.2f" jobname bi bf bstep
 
-    jobid = submitJob(env, prog, args, jobname, Run = true)   
+    jobid = submitJob(env, prog, args, jobname, 
+                            partitian = machine,
+                            processor = processor,
+                            gpu_memory = gpu_memory,
+                            cpu_per_task = cpu_per_task,
+                            Run = true, 
+                            Wait = Wait)
 end
\ No newline at end of file

jobs/jobsCMPO_XXZ_2D_helical.jl

@@ -5,35 +5,29 @@ phys_model = "XXZ_2D_helical"
 env = "/home/sliang/JuliaCode/JuliaCMPO"
 prog = env * "/jobs/CMPO_$(phys_model).jl"
 
-device = GPU
-machine = a100
+processor = GPU
+machine = p100
 gpu_memory = 80
-logtag = Dates.format(now(), "yyyy-mm-dd")*"-"*"$(device)"
+logtag = Dates.format(now(), "yyyy-mm-dd")*"-"*"$(processor)"
 
 Wait = nothing
-cpu_per_task = 8
+cpu_per_task = 1
 
-#tag = "2022-06-21"*"-"*device
-tag = Dates.format(now(), "yyyy-mm-dd")*"-"*"$(device)"
-beta1 = [0.1, 0.5]
-beta2 = [i for i in range(1.0, 5.0, step=1.0)]
-beta3 = [i for i in range(1.1, 1.9, step=0.2)]
-beta4 = [i for i in range(1.2, 1.9, step=0.2)]
+#tag = "2022-06-21"*"-"*processor
+tag = Dates.format(now(), "yyyy-mm-dd")*"-"*"$(processor)"
 
-#βlist = [i for i in range(1.1, 1.5, step=0.1)]
-#βlist = vcat(beta2, beta3, beta1)
-βlist = [5.0] 
+βlist = [0.5] 
 Jzlist = [1.0]
 Jxylist = [1.0]
-bondDlist = [64]
-widlist = [5]
+bondDlist = [32]
+widlist = [3]
 Continue = 0  #Continue > max_pow_step,  Continue = true
-max_pow_step = 1
+max_pow_step = 10
 
 
 
 #CREAT LOG FOLDER
-logdir = "/data/sliang/log/JuliaCMPO"
+logdir = "/data/sliang/JuliaCMPO/log"
 isdir(logdir) || mkdir(logdir)
 
 for bondD in bondDlist
@@ -46,7 +40,7 @@ for bondD in bondDlist
                     "max_pow_step"=>max_pow_step,
                     "Continue"=>Continue,
                     "tag"=>tag,
-                    "device"=> device
+                    "processor"=> Int(processor)
                     )
         jobname = logdir * "/" * phys_model
         isdir(jobname) || mkdir(jobname)
@@ -54,11 +48,12 @@ for bondD in bondDlist
         isdir(jobname) || mkdir(jobname)
         jobname = @sprintf "%s/bondD_%02i_%s" jobname bondD logtag
         isdir(jobname) || mkdir(jobname)
-        jobname = @sprintf "%s/beta_%.2f" jobname β
+        timetag = Dates.format(now(), "HH-MM-SS")
+        jobname = @sprintf "%s/beta_%.2f_%s" jobname β timetag
 
         jobid = submitJob(env, prog, args, jobname, 
                             partitian = machine,
-                            device = device,
+                            processor = processor,
                             gpu_memory = gpu_memory,
                             cpu_per_task = cpu_per_task,
                             Run = true, 

src/CMPSCompress.jl

@@ -2,7 +2,7 @@
     `compress_cmps`:
 """
 function compress_cmps(ψ0::AbstractCMPS{T, S, U}, χ::Integer, β::Real; 
-    options::CompressOptions=CompressOptions()) where {T,S,U}
+    options::CompressOptions=CompressOptions(trace_estimator=nothing)) where {T,S,U}
     @unpack (init, show_trace, mera_update_options, 
              optim_options, trace_estimator, processor) = options
     if show_trace 
@@ -33,7 +33,6 @@ function compress_cmps(ψ0::AbstractCMPS{T, S, U}, χ::Integer, β::Real;
             dQ = convert(Vector, diag(ψ.Q))
             R =  convert(Array, ψ.R)
 
-            T <: CuCMPS ? solver = gpu_solver : solver = cpu_solver
             loss() = -logfidelity(solver(CMPS_generate, diagm(dQ), R), ψ0, β, trace_estimator)
             p0, f, g! = optim_functions(loss, Params([dQ, R]))
 

src/CMPSOperations.jl

@@ -14,9 +14,9 @@ end
     Logarithm of the overlap of two CMPS:
         `log_overlap = ln(⟨ψl|ψr⟩)`
 """
-function log_overlap(ψl::AbstractCMPS, ψr::AbstractCMPS, β::Real, estimator::EstimatorType)
+function log_overlap(ψl::AbstractCMPS, ψr::AbstractCMPS, β::Real, trace_estimator::EstimatorType)
     K = CMPSMatrix(ψl, ψr)
-    return logtrexp(-β, K, estimator)
+    return logtrexp(-β, K, trace_estimator)
 end
 #log_overlap(ψl::AbstractCMPS, ψr::AbstractCMPS, β::Real) = log_overlap(ψl, ψr, β, nothing)
 
@@ -24,8 +24,8 @@ end
 """
     Norm of a CMPS
 """
-function norm(s::AbstractCMPS, β::Real, estimator::EstimatorType)
-    λ = log_overlap(s, s, β, estimator) 
+function norm(s::AbstractCMPS, β::Real, trace_estimator::EstimatorType)
+    λ = log_overlap(s, s, β, trace_estimator) 
     return exp(λ/2)
 end
 #norm(s::AbstractCMPS, β::Real) = norm(s, β, nothing)
@@ -34,8 +34,8 @@ end
 """
     Normalize a CMPS, i.e. ⟨ψ|ψ⟩ = 1
 """
-function normalize(s::AbstractCMPS{T, S, U}, β::Real, estimator::EstimatorType) where {T,S,U}
-    λ = log_overlap(s, s, β, estimator)/β
+function normalize(s::AbstractCMPS{T, S, U}, β::Real, trace_estimator::EstimatorType) where {T,S,U}
+    λ = log_overlap(s, s, β, trace_estimator)/β
     eye = λ/2 * Matrix{Float64}(I,size(s.Q))
     Q = s.Q - convert(S, eye)
     return CMPS_generate(Q, s.R)
@@ -48,19 +48,19 @@ end
         fidelity = ⟨ψ|T|r⟩/√(⟨ψ|ψ⟩)
         logfidelity(ψ, ψ0) = ln(Fd)
 """
-function logfidelity(ψ::T, ψ0::T, β::Real, esitmator::EstimatorType) where T<:AbstractCMPS
-    return log_overlap(ψ, ψ0, β, estimator) - 0.5*log_overlap(ψ, ψ, β, estimator)
+function logfidelity(ψ::T, ψ0::T, β::Real, trace_estimator::EstimatorType) where T<:AbstractCMPS
+    return log_overlap(ψ, ψ0, β, trace_estimator) - 0.5*log_overlap(ψ, ψ, β, trace_estimator)
 end
 #logfidelity(ψ::T, ψ0::T, β::Real) where T<:AbstractCMPS = logfidelity(ψ, ψ0, β, nothing)
 
 
-function fidelity(ψ::T, ψ0::T, β::Real, estimator::EstimatorType; 
+function fidelity(ψ::T, ψ0::T, β::Real, trace_estimator::EstimatorType; 
                   Normalize::Bool = false) where T<:AbstractCMPS
     if Normalize
-        ψ = normalize(ψ, β, estimator)
-        ψ0 = normalize(ψ0, β, estimator)
+        ψ = normalize(ψ, β, trace_estimator)
+        ψ0 = normalize(ψ0, β, trace_estimator)
     end
-    return logfidelity(ψ, ψ0, β, estimator) |> exp
+    return logfidelity(ψ, ψ0, β, trace_estimator) |> exp
 end
 #fidelity(ψ::T, ψ0::T, β::Real; Normalize::Bool = false) where T<:AbstractCMPS =
 #    fidelity(ψ, ψ0, β, nothing; Normalize = Normalize)

src/JuliaCMPO.jl

@@ -37,6 +37,7 @@ export MeraUpdateOptions,
        MeraUpdateResult
 #structs/CMPSCompress.jl
 export CompressOptions, CompressResult
+#export update_processor
 
 #utilities
 export ⊗
@@ -109,7 +110,8 @@ include("./PhysicalQuantities/Correlations.jl")
 include("./structs/EvaluateOptions.jl")
 include("evaluate.jl")
 
-
+include("./rrule/ConstructorAdjoint.jl")
+include("./rrule/accum.jl")
 include("./rrule/logtrexp.jl")
 
 

src/MeraUpdate.jl

@@ -14,7 +14,7 @@ end
 
 
 function adaptive_mera_update(ψ0::AbstractCMPS, χ::Integer, β::Real; 
-    options::MeraUpdateOptions = MeraUpdateOptions())
+    options::MeraUpdateOptions = MeraUpdateOptions(trace_estimator=nothing))
     @unpack (atol, ldiff_tol, maxiter, interpolate,
             store_trace, show_trace, trace_estimator) = options
     step = 1
@@ -22,7 +22,8 @@ function adaptive_mera_update(ψ0::AbstractCMPS, χ::Integer, β::Real;
     logfidelity0 = 9.9e9
     loss(p_matrix) = logfidelity(project(ψ0, p_matrix), ψ0, β, trace_estimator)
 
-    _, v = eigensolver(ψ0.Q |> symmetrize)
+    Q = symmetrize(ψ0.Q)
+    _, v = eigensolver(Q)
     p_current = v[:, end-χ+1:end]
     
     loss_previous = 9.9e9

src/OptimFunctions.jl

@@ -7,6 +7,7 @@ veclength(x) = length(x)
 Base.zeros(grads::Zygote.Grads) = zeros(veclength(grads))
 Base.zeros(pars::Zygote.Params) = zeros(veclength(pars))
 
+
 """
     optim_function(loss, pars): Return two functions(loss function, gradient function)
     and p0, a vectorized version of pars.

src/PhysicalQuantities/Correlations.jl

@@ -5,7 +5,8 @@
 function correlation_2time(τ::Number, A::AbstractMatrix,B::AbstractMatrix,
                             ψl::AbstractCMPS, ψr::AbstractCMPS, W::AbstractCMPO, β::Real)
     K = ψl * W * ψr 
-    e, v = symeigen(K)
+    K = symmetrize(K)
+    e, v = eigensolver(K)
     min = minimum(e); e = e .- min
     A = v' * A * v
     B = v' * B * v

src/PhysicalQuantities/PhysicalModels.jl

 #module PhysicalModels
 """
-    Data structure of a physical model
+    Data structure of a physical model 
     PhysModel:
         `Tmatrix`: local transfer matrix
         `phy_dim`: bond dimension of physical legs
@@ -9,21 +9,22 @@
               Note that in general, `Ut` not necessarily be unitary, but be invertible, 
               however, in the limiting cases we know right now, they are unitary.
 """
-struct PhysModel{T<:AbstractCMPO}
-    Tmatrix::T 
+@with_kw struct PhysModel{Tm<:AbstractCMPO, 
+                          Tu<:Union{AbstractMatrix, Nothing}}
+    Tmatrix::Tm 
     phy_dim::Int64  
     vir_dim::Int64  
-    Ut::Union{AbstractMatrix, Nothing} 
+    Ut::Tu
 end
 
 
 """
     Ising type CMPO: H = -J ôl ôr block
 """
-function Ising_CMPO(J::Real, ol::AbstractArray, or::AbstractArray, wid::Integer = 1)
+function Ising_CMPO(J::Real, ol::AbstractArray{T}, or::AbstractArray{T}, wid::Integer = 1) where {T}
     sgn = sign(J); val =√(abs(J))
     o2 = zeros(2, 2)
-    i2 = Matrix{Float64}(I, 2, 2)
+    i2 = Matrix{T}(I, 2, 2)
     #construct L and R
     ol = val * ol; or = sgn * val * or
     L = ol; R = or
@@ -42,7 +43,7 @@ function Ising_CMPO(J::Real, ol::AbstractArray, or::AbstractArray, wid::Integer
         i == 1 ? P = pcol : P = cat(P, pcol, dims = 4)
     end
         
-    return CMPO(zeros(2,2), R, L, P)
+    return CMPO(zeros(T,2,2), R, L, P)
 end
 
 
@@ -61,18 +62,18 @@ end
 """
     Expand cmpo
 """
-function expand_cmpo(o::CMPO)
+function expand_cmpo(o::AbstractCMPO{T,S,U,V}) where {T,S,U,V}
     R = √(0.5) * cat(o.R, o.L, dims = 3)
     L = √(0.5) * cat(o.L, o.R, dims = 3)
 
-    if length(size(o.P)) == 2
-        oP = reshape(o.P, size(o.P)[1], size(o.P)[2], 1, 1)
+    if V <: AbstractMatrix
+        oP = reshape(o.P, size(o.P,1), size(o.P,2), 1, 1)
     else
         oP = o.P
     end
 
     pl = cat(oP, zeros(size(oP)), dims = 3)
-    pr = cat(zeros(size(oP)), ein"ijkl->ijlk"(oP), dims = 3)
+    pr = cat(zeros(T,size(oP)), ein"ijkl->ijlk"(oP), dims = 3)
     P = cat(pl, pr, dims = 4)
     
     return CMPO(o.Q, R, L, P)
@@ -82,17 +83,19 @@ end
 """
     cat CMPO blocks
 """
-function cat(o1::CMPO, o2::CMPO)
-    pd = size(o1.Q)[1]
-    length(size(o1.P)) == 2 ? vd1 = 1 : vd1 = size(o1.P)[3]
-    length(size(o2.P)) == 2 ? vd2 = 1 : vd2 = size(o2.P)[3]
+function cat(o1::AbstractCMPO{T,S1,U1,V1}, 
+             o2::AbstractCMPO{T,S2,U2,V2}
+             ) where {T,S1,U1,V1, S2,U2,V2}
+    pd = size(o1.Q,1)
+    V1<:AbstractMatrix ? vd1 = 1 : vd1 = size(o1.P,3)
+    V2<:AbstractMatrix ? vd2 = 1 : vd2 = size(o2.P,3)
 
-    Q = zeros(pd, pd)
+    Q = zeros(T, pd, pd)
     R = cat(o1.R, o2.R, dims = 3)
     L = cat(o1.L, o2.L, dims = 3)
 
-    pl = cat(o1.P, zeros(pd, pd, vd2, vd1), dims = 3)
-    pr = cat(zeros(pd, pd, vd1, vd2), o2.P, dims = 3)
+    pl = cat(o1.P, zeros(T, pd, pd, vd2, vd1), dims = 3)
+    pr = cat(zeros(T, pd, pd, vd1, vd2), o2.P, dims = 3)
     P = cat(pl, pr, dims = 4)
 
     return CMPO(Q, R, L, P)
@@ -128,7 +131,7 @@ function XYmodel()
     Q = zeros(2, 2)
     P = zeros(2, 2, 2, 2)
     Tmatrix = CMPO(Q,R,L,P)
-    Ut = Matrix(1.0I, 2, 2)
+    Ut = Matrix{Float64}(I, 2, 2)
     return PhysModel(Tmatrix, 2, 2, Ut)
 end
 
@@ -171,7 +174,7 @@ function XXZmodel(Δ::Real)
         Q = zeros(2, 2)
         P = zeros(2, 2, 3, 3)
         Tmatrix = CMPO(Q,R,L,P)
-        Ut = Matrix(1.0I, 3, 3) 
+        Ut = Matrix{Float64}(I, 3, 3) 
         return PhysModel(Tmatrix, 2, 3, Ut)
     end
 end

src/PhysicalQuantities/Thermaldynamic.jl

@@ -2,12 +2,12 @@
 """
     make operator ⊢o⊣
 """
-function make_operator(Op::AbstractMatrix{T}, dim::Int64) where T
+function make_operator(Op::AbstractMatrix{T}, dim::Int64) where {T}
     eye = Matrix{T}(I, dim, dim)
     return eye ⊗ Op ⊗ eye
 end
 
-function make_operator(Op::AbstractMatrix{T}, ψ::AbstractCMPS) where T
+function make_operator(Op::AbstractMatrix{T}, ψ::AbstractCMPS{T,S,U}) where {T,S,U}
     eye = Matrix{T}(I, size(ψ.Q))
     return eye ⊗ Op ⊗ eye
 end
@@ -17,7 +17,8 @@ end
     The thermal average of local opeartors ⊢o⊣ with respect to K = ψl * W * ψr
 """
 function thermal_average(Op::AbstractMatrix, ψl::AbstractCMPS, ψr::AbstractCMPS, W::AbstractCMPO, β::Real)
-    K = ψl * W * ψr 
+    K = ψl * W * ψr |> Matrix |> symmetrize
+
     e, v = eigensolver(K)
     e = -β*e
     m = maximum(e); e = e .- m
@@ -34,7 +35,8 @@ thermal_average(Op::AbstractMatrix, ψ::AbstractCMPS, W::AbstractCMPO, β::Real)
     The thermal average of local opeartors ⊢o⊣ with respect to K = ψ * ψ
 """
 function thermal_average(Op::AbstractMatrix, ψl::AbstractCMPS, ψr::AbstractCMPS, β::Real)
-    K = ψl * ψr 
+    K = ψl * ψr |> Matrix |> symmetrize
+    
     e, v = eigensolver(K)
     e = -β*e
     m = maximum(e) ; e = e .- m
@@ -43,25 +45,27 @@ function thermal_average(Op::AbstractMatrix, ψl::AbstractCMPS, ψr::AbstractCMP
     num = exp.(e) .* diag(Op) |> sum
     return num/den
 end
-thermal_average(Op::AbstractMatrix, ψ::AbstractCMPS, β::Real) =thermal_average(Op, ψ, ψ, β)
+thermal_average(Op::AbstractMatrix, ψ::AbstractCMPS, β::Real) = thermal_average(Op, ψ, ψ, β)
 
 
 """
     Free energy: F = -1/(βL) lnZ
 """
-function free_energy(ψl::AbstractCMPS, ψr::AbstractCMPS, W::AbstractCMPO, β::Real)
-    res = log_overlap(ψl, W * ψr, β) - log_overlap(ψl, ψr, β)
+function free_energy(ψl::AbstractCMPS, ψr::AbstractCMPS, 
+                     W::AbstractCMPO, β::Real, trace_estimator::EstimatorType = nothing)
+    res = log_overlap(ψl, W * ψr, β, trace_estimator) - log_overlap(ψl, ψr, β, trace_estimator)
     return -res/β
 end
-free_energy(ψ::AbstractCMPS, W::AbstractCMPO, β::Real) = free_energy(ψ, ψ, W, β)
+free_energy(ψ::AbstractCMPS, W::AbstractCMPO, β::Real, trace_estimator::EstimatorType = nothing) = 
+    free_energy(ψ, ψ, W, β, trace_estimator)
 
 
 """
     Energy density: E = -∂lnZ/∂β
 """
 function energy(ψl::AbstractCMPS, ψr::AbstractCMPS, W::AbstractCMPO, β::Real)
-    K = ψl * W * ψr 
-    H = ψl * ψr 
+    K = ψl * W * ψr |> Matrix |> symmetrize
+    H = ψl * ψr |> Matrix |> symmetrize
     res = thermal_average(K, ψl, ψr, W, β) - thermal_average(H, ψl, ψr, β)
     return res
 end
@@ -74,8 +78,8 @@ energy(ψ::AbstractCMPS, W::AbstractCMPO, β::Real) = energy(ψ, ψ, W, β)
 function specific_heat(ψl::AbstractCMPS, ψr::AbstractCMPS, W::AbstractCMPO, β::Real; 
                         method::Symbol = :adiff)
     if method == :adiff
-        K = ψl * W * ψr 
-        H = ψl * ψr 
+        K = ψl * W * ψr |> Matrix |> symmetrize
+        H = ψl * ψr |> Matrix |> symmetrize
         K2 = K * K
         H2 = H * H
         c = thermal_average(K2, ψl, ψr, W, β) - thermal_average(K, ψl, ψr, W, β)^2

src/eigensolver.jl

 @reexport import KrylovKit.eigsolve
 import LinearAlgebra.Eigen
 
-function eigensolver(M::CMPSMatrix{Ts,T,S,U}; 
-                     which::Symbol=:SR, 
-                     ishermitian = true, 
-                     kwargs...) where {Ts,T,S,U} 
-    N = size(M,1)
-    Ts <: CMPS ? x0 = rand(T,N) : x0 = CUDA.rand(T, N)
-    vals, vecs, _ = eigsolve(M, x0, N, which, ishermitian = ishermitian, krylovdim = N, kwargs...)
-    vals = convert(typeof(x0), vals)
-    vecs = hcat(vecs[1:N]...)
-    return Eigen(vals, vecs)
+function eigensolver(M::CMPSMatrix{Ts,T,S,U}) where {Ts,T,S,U} 
+    res = Matrix(M) |> symmetrize
+    return eigensolver(res)
 end
 
-#function eigsolve(M::CMPSMatrix{Ts,T,S,U}, x0, howmany, which; args...) where {Ts,T,S,U}
-#    Ts <: CMPS ? x0 = rand(size(M,1)) : x0 = CUDA.rand(size(M,1))
-#    return eigsolve(M, x0, howmany, which, T, ishermitian = true, args...)
-#end
-
-    
\ No newline at end of file
+#end
\ No newline at end of file

src/evaluate.jl

@@ -19,22 +19,18 @@ end
     Evaluate PhysModel m when its transfer matrix is hermitian
 """
 function hermitian_evaluate(m::PhysModel, bondD::Integer, β::Real, ResultFolder::String; 
-            options::EvaluateOptions=EvaluateOptions())
+            options::EvaluateOptions=EvaluateOptions(trace_estimator=nothing))
     @unpack (init, processor, trace_estimator, 
-            compress_options, optim_options, tag) = options
-    if trace_estimator !== nothing
-        new_options = FTLMOptions(trace_estimator.options, processor=processor)
-        trace_estimator = TraceEstimator(trace_estimator.estimator, new_options)
-        compress_options = CompressOptions(compress_options, trace_estimator = trace_estimator)
-    end
+            compress_options, optim_options, tag, show_trace) = options
     solver = solver_function(processor)
 
     """
     ResultFolder: classified by Model parameters(interaction, width), store CMPS and Obsv files
     CMPSResultFolder: CMPS information, classified by bond dimension
     """
-    CMPSResultFolder = @sprintf "%s/bondD_%02i_CMPS_%s" ResultFolder bondD tag
-    OptResultFolder = @sprintf "%s/bondD_%02i_Opt_%s" ResultFolder bondD tag
+    trace_estimator === nothing ? estimator_name = "nothing" : estimator_name = string(trace_estimator.estimator)
+    CMPSResultFolder = @sprintf "%s/bondD_%02i_CMPS_%s_%s" ResultFolder bondD estimator_name tag
+    OptResultFolder = @sprintf "%s/bondD_%02i_Opt_%s_%s" ResultFolder bondD estimator_name tag
     isdir(ResultFolder) || mkdir(ResultFolder)
     isdir(CMPSResultFolder) || mkdir(CMPSResultFolder) 
     isdir(OptResultFolder) || mkdir(OptResultFolder)
@@ -49,10 +45,10 @@ function hermitian_evaluate(m::PhysModel, bondD::Integer, β::Real, ResultFolder
     dQ = convert(Vector, diag(ψ.Q))
     R  = convert(Array, ψ.R)
 
-    loss() = free_energy(solver(CMPS_generate, diagm(dQ), R), Tm, β)
+    loss() = free_energy(solver(CMPS_generate, diagm(dQ), R), Tm, β, trace_estimator)
     F_initial = loss()
 
-    p0, f, g! = optim_functions(loss, Params([dQ, R]))
+    p0, f, g! = optim_functions(loss, Zygote.Params([dQ, R]))
     optim_result = Optim.optimize(f, g!, p0, LBFGS(), optim_options)
     F_final = minimum(optim_result)
 
@@ -86,14 +82,9 @@ end
 """
 function non_hermitian_evaluate(m::PhysModel, bondD::Integer, β::Real, ResultFolder::String; 
                                 options::EvaluateOptions=EvaluateOptions())
-    @unpack (init, solver, trace_estimator, 
+    @unpack (init, processor, trace_estimator, 
             compress_options, optim_options, 
-            Continue, max_pow_step, group, tag) = options
-    if trace_estimator !== nothing
-        new_options = FTLMOptions(trace_estimator.options, processor=processor)
-        trace_estimator = TraceEstimator(trace_estimator.estimator, new_options)
-        compress_options = CompressOptions(compress_options, trace_estimator = trace_estimator)
-    end
+            Continue, max_pow_step, group, tag, show_trace) = options
     solver = solver_function(processor)
     
     """
@@ -108,7 +99,8 @@ function non_hermitian_evaluate(m::PhysModel, bondD::Integer, β::Real, ResultFo
         g+=1
     end
 
-    CMPSResultFolder = @sprintf "%s/bondD_%02i_CMPS_%s" ResultFolder bondD tag
+    trace_estimator === nothing ? estimator_name = "nothing" : estimator_name = string(trace_estimator.estimator)
+    CMPSResultFolder = @sprintf "%s/bondD_%02i_CMPS_%s_%s" ResultFolder bondD estimator_name tag
     ChkpFolder = @sprintf "%s/CheckPoint_beta_%.2f" CMPSResultFolder β
     ChkpPsiFolder = @sprintf "%s/cmps_psi" ChkpFolder
     ChkpLpsiFolder = @sprintf "%s/cmps_Lpsi" ChkpFolder

src/logtrexp.jl

@@ -17,7 +17,7 @@ function logtrexp(t::Real, M, trace_estimator::TraceEstimator)
     @unpack processor = options
     processor == CPU ? x0 = rand(size(M,1)) : x0 = CUDA.rand(Float64, size(M,1))
     e0, _, _ = eigsolve(M, x0, 1, which, ishermitian = true)
-    e1 = e0[1]
+    e1 = e0[1] #typeof(e0) = Vector{Float64}
     expr = e -> exp(t * (e - e1))
 
     options = FTLMOptions(options, which = which)