Merge branch 'master' of github.com:quantumghent/PEPSKit.jl

Author: pbrehmer

examples/heisenberg_evol/simpleupdate.jl

@@ -7,9 +7,7 @@ trscheme_env = truncerr(1e-10) & truncdim(χenv)
 Nr, Nc = 2, 2
 # Heisenberg model Hamiltonian
 # (already only includes nearest neighbor terms)
-ham = heisenberg_XYZ(ComplexF64, symm, InfiniteSquare(Nr, Nc); Jx=1.0, Jy=1.0, Jz=1.0)
-# convert to real tensors
-ham = LocalOperator(ham.lattice, Tuple(ind => real(op) for (ind, op) in ham.terms)...)
+ham = real(heisenberg_XYZ(ComplexF64, symm, InfiniteSquare(Nr, Nc); Jx=1.0, Jy=1.0, Jz=1.0))
 
 # random initialization of 2x2 iPEPS with weights and CTMRGEnv (using real numbers)
 if symm == Trivial

src/algorithms/optimization/peps_optimization.jl

@@ -112,9 +112,10 @@ The optimization parameters can be supplied via the keyword arguments or directl
 * `tol::Real=$(Defaults.optimizer_tol)` : Overall tolerance for gradient norm convergence of the optimizer. Sets related tolerance such as the boundary and boundary-gradient tolerances to sensible defaults unless they are explictly specified.
 * `verbosity::Int=1` : Overall output information verbosity level, should be one of the following:
     0. Suppress all output
-    1. Optimizer output and warnings
-    2. Additionally print boundary information
-    3. All information including AD debug outputs
+    1. Only print warnings
+    2. Initialization and convergence info
+    3. Iteration info
+    4. Debug info including AD outputs
 * `reuse_env::Bool=$(Defaults.reuse_env)` : If `true`, the current optimization step is initialized on the previous environment, otherwise a random environment is used.
 * `symmetrization::Union{Nothing,SymmetrizationStyle}=nothing` : Accepts `nothing` or a `SymmetrizationStyle`, in which case the PEPS and PEPS gradient are symmetrized after each optimization iteration.
 * `(finalize!)=OptimKit._finalize!` : Inserts a `finalize!` function call after each optimization step by utilizing the `finalize!` kwarg of `OptimKit.optimize`. The function maps `(peps, env), f, g = finalize!((peps, env), f, g, numiter)`.

src/algorithms/select_algorithm.jl

@@ -21,29 +21,29 @@ function select_algorithm(
     ::typeof(fixedpoint),
     env₀::CTMRGEnv;
     tol=Defaults.optimizer_tol, # top-level tolerance
-    verbosity=2, # top-level verbosity
+    verbosity=3, # top-level verbosity
     boundary_alg=(;),
     gradient_alg=(;),
     optimizer_alg=(;),
     kwargs...,
 )
     # adjust CTMRG tols and verbosity
     if boundary_alg isa NamedTuple
-        defaults = (; verbosity=verbosity ≤ 1 ? -1 : verbosity, tol=1e-4tol)
+        defaults = (; verbosity=verbosity ≤ 3 ? -1 : 3, tol=1e-4tol)
         boundary_kwargs = merge(defaults, boundary_alg)
         boundary_alg = select_algorithm(leading_boundary, env₀; boundary_kwargs...)
     end
 
     # adjust gradient verbosity
     if gradient_alg isa NamedTuple
         # TODO: check this:
-        defaults = (; verbosity=verbosity ≤ 2 ? -1 : 3, tol=1e-2tol)
+        defaults = (; verbosity=verbosity ≤ 3 ? -1 : 3, tol=1e-2tol)
         gradient_alg = merge(defaults, gradient_alg)
     end
 
     # adjust optimizer tol and verbosity
     if optimizer_alg isa NamedTuple
-        defaults = (; tol, verbosity=verbosity ≤ 1 ? -1 : 3)
+        defaults = (; tol, verbosity)
         optimizer_alg = merge(defaults, optimizer_alg)
     end
 

src/operators/localoperator.jl

@@ -14,7 +14,7 @@ and the terms are stored as a tuple of pairs of indices and operators.
 # Constructors
 
     LocalOperator(lattice::Matrix{S}, terms::Pair...)
-    LocalOperator{T,S}(lattice::Matrix{S}, terms::T) where {T,S} # expert mode
+    LocalOperator{T,S}(lattice::Matrix{S}, terms::T) where {T,S}
 
 # Examples
 
@@ -31,6 +31,7 @@ struct LocalOperator{T<:Tuple,S}
         # Check if the indices of the operator are valid with themselves and the lattice
         for (inds, operator) in terms
             @assert operator isa AbstractTensorMap
+            @assert eltype(inds) <: CartesianIndex
             @assert numout(operator) == numin(operator) == length(inds)
             @assert spacetype(operator) == S
 
@@ -52,7 +53,12 @@ function LocalOperator(
     relevant_terms = []
     for inds in unique(allinds)
         operator = sum(alloperators[findall(==(inds), allinds)])
-        norm(operator) > atol && push!(relevant_terms, inds => operator)
+        cinds = if !(eltype(inds) <: CartesianIndex) # force indices to be CartesianIndices
+            map(CartesianIndex, inds)
+        else
+            inds
+        end
+        norm(operator) > atol && push!(relevant_terms, cinds => operator)
     end
 
     terms_tuple = Tuple(relevant_terms)
@@ -90,6 +96,15 @@ function Base.repeat(O::LocalOperator, m::Int, n::Int)
     return LocalOperator(lattice, terms...)
 end
 
+# Real and imaginary part
+# -----------------------
+function Base.real(O::LocalOperator)
+    return LocalOperator(O.lattice, (ind => real(op) for (ind, op) in O.terms)...)
+end
+function Base.imag(O::LocalOperator)
+    return LocalOperator(O.lattice, (ind => imag(op) for (ind, op) in O.terms)...)
+end
+
 # Linear Algebra
 # --------------
 function Base.:*(α::Number, O::LocalOperator)

src/states/infiniteweightpeps.jl

@@ -267,6 +267,21 @@ function InfinitePEPS(peps::InfiniteWeightPEPS)
     )
 end
 
+"""
+    InfiniteWeightPEPS(peps::InfinitePEPS)
+
+Create `InfiniteWeightPEPS` from `InfinitePEPS` by initializing the bond weights as identity matrices of element type `Float64`.
+"""
+function InfiniteWeightPEPS(peps::InfinitePEPS)
+    Nr, Nc = size(peps)
+    weights = map(Iterators.product(1:2, 1:Nr, 1:Nc)) do (d, r, c)
+        V = (d == 1 ? domain(peps[r, c])[2] : domain(peps[r, c])[1])
+        @assert !isdual(V)
+        DiagonalTensorMap(ones(reduceddim(V)), V)
+    end
+    return InfiniteWeightPEPS(peps.A, SUWeight(weights))
+end
+
 """
     mirror_antidiag(peps::InfiniteWeightPEPS)
 

test/examples/heisenberg.jl

@@ -59,10 +59,11 @@ end
     for ind in CartesianIndices(wpeps.vertices)
         wpeps.vertices[ind] /= norm(wpeps.vertices[ind], Inf)
     end
-    # Heisenberg model Hamiltonian (already only includes nearest neighbor terms)
+    # Heisenberg model Hamiltonian
     ham = heisenberg_XYZ(InfiniteSquare(N1, N2); Jx=1.0, Jy=1.0, Jz=1.0)
-    # convert to real tensors
-    ham = LocalOperator(ham.lattice, Tuple(ind => real(op) for (ind, op) in ham.terms)...)
+    # assert imaginary part is zero
+    @assert length(imag(ham).terms) == 0
+    ham = real(ham)
 
     # simple update
     dts = [1e-2, 1e-3, 1e-3, 1e-4]