From ca7418d2cd0e460caaa58fe736e6a8464ef54205 Mon Sep 17 00:00:00 2001 From: mtfishman Date: Sat, 11 Jan 2025 19:51:59 -0500 Subject: [PATCH] Work towards getting non-QN DMRG working --- examples/dmrg/Project.toml | 2 + src/ITensorMPS.jl | 1 - src/abstractmps.jl | 5 +- src/abstractprojmpo/projmposum.jl | 4 +- src/mpo.jl | 2 +- src/opsum_to_mpo/opsum_to_mpo.jl | 92 ++++++++++++++++++++++++ src/opsum_to_mpo/opsum_to_mpo_generic.jl | 30 ++++++-- src/solvers/timedependentsum.jl | 5 +- 8 files changed, 129 insertions(+), 12 deletions(-) diff --git a/examples/dmrg/Project.toml b/examples/dmrg/Project.toml index e6b9177..7b391c3 100644 --- a/examples/dmrg/Project.toml +++ b/examples/dmrg/Project.toml @@ -1,6 +1,8 @@ [deps] +ITensorBase = "4795dd04-0d67-49bb-8f44-b89c448a1dc7" ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +NamedDimsArrays = "60cbd0c0-df58-4cb7-918c-6f5607b73fde" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" diff --git a/src/ITensorMPS.jl b/src/ITensorMPS.jl index 30add93..1c42b83 100644 --- a/src/ITensorMPS.jl +++ b/src/ITensorMPS.jl @@ -43,5 +43,4 @@ include("solvers/reducedlinearproblem.jl") include("solvers/linsolve.jl") include("solvers/expand.jl") include("lib/Experimental/src/Experimental.jl") -## include("lib/ITensorMPSNamedDimsArraysExt/src/ITensorMPSNamedDimsArraysExt.jl") end diff --git a/src/abstractmps.jl b/src/abstractmps.jl index 9e823c6..014bdd1 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -2,7 +2,8 @@ using ITensors: @Algorithm_str using IsApprox: Approx, IsApprox -using ITensors: ITensors, ITensor, Index, commonind, commoninds, uniqueind, uniqueinds +using ITensors: + ITensors, ITensor, Index, commonind, commoninds, dag, hasqns, uniqueind, uniqueinds ## using NDTensors: NDTensors, using_auto_fermion, scalartype, tensor using ITensors.Ops: Prod ## using ITensors.QuantumNumbers: QuantumNumbers, removeqn @@ -2293,7 +2294,7 @@ end Return true if the MPS or MPO has tensors which carry quantum numbers. """ -hasqns(M::AbstractMPS) = any(hasqns, data(M)) +ITensors.hasqns(M::AbstractMPS) = any(hasqns, data(M)) # Trait type version of hasqns # Note this is not inferrable, so hasqns would be preferred diff --git a/src/abstractprojmpo/projmposum.jl b/src/abstractprojmpo/projmposum.jl index 3d5d5d9..571c44d 100644 --- a/src/abstractprojmpo/projmposum.jl +++ b/src/abstractprojmpo/projmposum.jl @@ -1,8 +1,8 @@ -using Compat: allequal +using ITensors.Ops: Ops abstract type AbstractSum end -terms(sum::AbstractSum) = sum.terms +Ops.terms(sum::AbstractSum) = sum.terms function set_terms(sum::AbstractSum, terms) return error("Please implement `set_terms` for the `AbstractSum` type `$(typeof(sum))`.") diff --git a/src/mpo.jl b/src/mpo.jl index 4e19724..828ee0b 100644 --- a/src/mpo.jl +++ b/src/mpo.jl @@ -1,7 +1,7 @@ using Adapt: adapt using LinearAlgebra: dot using Random: Random, AbstractRNG -using ITensors: @ts_str, Algorithm, outer +using ITensors: @ts_str, Algorithm, dag, outer using ITensors.Ops: OpSum using ITensors.SiteTypes: SiteTypes, siteind, siteinds diff --git a/src/opsum_to_mpo/opsum_to_mpo.jl b/src/opsum_to_mpo/opsum_to_mpo.jl index 9142e02..15e855c 100644 --- a/src/opsum_to_mpo/opsum_to_mpo.jl +++ b/src/opsum_to_mpo/opsum_to_mpo.jl @@ -1,4 +1,96 @@ ## using NDTensors: using_auto_fermion +using ITensors: dim + +replace_nothing(::Nothing, y) = y +replace_nothing(x, y) = x + +default_mindim(a) = true +default_use_absolute_cutoff(a) = false +default_use_relative_cutoff(a) = true + +function truncate!( + P::AbstractVector; + mindim=nothing, + maxdim=nothing, + cutoff=nothing, + use_absolute_cutoff=nothing, + use_relative_cutoff=nothing, +) + mindim = replace_nothing(mindim, default_mindim(P)) + maxdim = replace_nothing(maxdim, length(P)) + cutoff = replace_nothing(cutoff, typemin(eltype(P))) + use_absolute_cutoff = replace_nothing(use_absolute_cutoff, default_use_absolute_cutoff(P)) + use_relative_cutoff = replace_nothing(use_relative_cutoff, default_use_relative_cutoff(P)) + + origm = length(P) + docut = zero(eltype(P)) + + #if P[1] <= 0.0 + # P[1] = 0.0 + # resize!(P, 1) + # return 0.0, 0.0 + #end + + if origm == 1 + docut = abs(P[1]) / 2 + return zero(eltype(P)), docut + end + + s = sign(P[1]) + s < 0 && (P .*= s) + + #Zero out any negative weight + for n in origm:-1:1 + (P[n] >= zero(eltype(P))) && break + P[n] = zero(eltype(P)) + end + + n = origm + truncerr = zero(eltype(P)) + while n > maxdim + truncerr += P[n] + n -= 1 + end + + if use_absolute_cutoff + #Test if individual prob. weights fall below cutoff + #rather than using *sum* of discarded weights + while P[n] <= cutoff && n > mindim + truncerr += P[n] + n -= 1 + end + else + scale = one(eltype(P)) + if use_relative_cutoff + scale = sum(P) + (scale == zero(eltype(P))) && (scale = one(eltype(P))) + end + + #Continue truncating until *sum* of discarded probability + #weight reaches cutoff reached (or m==mindim) + while (truncerr + P[n] <= cutoff * scale) && (n > mindim) + truncerr += P[n] + n -= 1 + end + + truncerr /= scale + end + + if n < 1 + n = 1 + end + + if n < origm + docut = (P[n] + P[n + 1]) / 2 + if abs(P[n] - P[n + 1]) < eltype(P)(1e-3) * P[n] + docut += eltype(P)(1e-3) * P[n] + end + end + + s < 0 && (P .*= s) + resize!(P, n) + return truncerr, docut +end # `ValType::Type{<:Number}` is used instead of `ValType::Type` for efficiency, possibly due to increased method specialization. # See https://github.com/ITensor/ITensors.jl/pull/1183. diff --git a/src/opsum_to_mpo/opsum_to_mpo_generic.jl b/src/opsum_to_mpo/opsum_to_mpo_generic.jl index 009a699..1d9a0cb 100644 --- a/src/opsum_to_mpo/opsum_to_mpo_generic.jl +++ b/src/opsum_to_mpo/opsum_to_mpo_generic.jl @@ -1,7 +1,28 @@ ## using NDTensors: using_auto_fermion -using ITensors.Ops: Ops, Op, OpSum, Scaled, Sum, coefficient +using ITensors.Ops: + Ops, Op, OpSum, Scaled, Sum, argument, coefficient, site, sites, terms, which_op using ITensors.SiteTypes: has_fermion_string, op +""" + parity_sign(P) + +Given an array or tuple of integers representing +a permutation or a subset of a permutation, +compute the parity sign defined as -1 for a +permutation consisting of an odd number of swaps +and +1 for an even number of swaps. This +implementation uses an O(n^2) algorithm and is +intended for small permutations only. +""" +function parity_sign(P)::Int + L = length(P) + s = +1 + for i in 1:L, j in (i + 1):L + s *= sign(P[j] - P[i]) + end + return s +end + # TODO: Deprecate. const AutoMPO = OpSum @@ -167,7 +188,7 @@ function sorteachterm(os::OpSum, sites) os = copy(os) for (j, t) in enumerate(os) - if maximum(ITensors.sites(t)) > length(sites) + if maximum(Ops.sites(t)) > length(sites) error( "The OpSum contains a term $t that extends beyond the number of sites $(length(sites)).", ) @@ -188,7 +209,8 @@ function sorteachterm(os::OpSum, sites) t_parity = +1 for n in reverse(1:Nt) site_n = only(site(t[n])) - if !using_auto_fermion() && (t_parity == -1) && (site_n < prevsite) + ## if !using_auto_fermion() && (t_parity == -1) && (site_n < prevsite) + if (t_parity == -1) && (site_n < prevsite) # Insert local piece of Jordan-Wigner string emanating # from fermionic operators to the right # (Remaining F operators will be put in by svdMPO) @@ -212,7 +234,7 @@ function sorteachterm(os::OpSum, sites) filter!(!iszero, perm) # and account for anti-commuting, fermionic operators # during above sort; put resulting sign into coef - t *= ITensors.parity_sign(perm) + t *= parity_sign(perm) terms(os)[j] = t end diff --git a/src/solvers/timedependentsum.jl b/src/solvers/timedependentsum.jl index f2998de..d556a7d 100644 --- a/src/solvers/timedependentsum.jl +++ b/src/solvers/timedependentsum.jl @@ -1,4 +1,5 @@ using ITensors: ITensor, inds, permute +using ITensors.Ops: Ops # Represents a time-dependent sum of terms: # @@ -10,7 +11,7 @@ struct TimeDependentSum{Coefficients,Terms} end coefficients(expr::TimeDependentSum) = expr.coefficients -terms(expr::TimeDependentSum) = expr.terms +Ops.terms(expr::TimeDependentSum) = expr.terms function Base.copy(expr::TimeDependentSum) return TimeDependentSum(coefficients(expr), copy.(terms(expr))) end @@ -51,7 +52,7 @@ struct ScaledSum{Coefficients,Terms} end coefficients(expr::ScaledSum) = expr.coefficients -terms(expr::ScaledSum) = expr.terms +Ops.terms(expr::ScaledSum) = expr.terms # Apply the scaled sum of terms: #