Work towards getting non-QN DMRG working

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"

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

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

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))`.")

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
 

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.

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
 

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:
 #