Merge remote-tracking branch 'upstream/master' into bp

Author: Yue-Zhengyuan

.github/workflows/FormatCheck.yml

@@ -11,4 +11,6 @@ on:
 jobs:
   format-check:
     name: "Format Check"
-    uses: "QuantumKitHub/.github/.github/workflows/formatcheck.yml@main"
+    uses: "QuantumKitHub/.github/.github/workflows/formatcheck.yml@main"
+    with: 
+      juliaformatter-version: "2"

docs/make.jl

@@ -41,14 +41,15 @@ mathengine = MathJax3(
 )
 
 # examples pages
-examples_optimization =
-    joinpath.(["heisenberg", "bose_hubbard", "xxz", "fermi_hubbard"], Ref("index.md"))
-examples_time_evolution =
-    joinpath.(["heisenberg_su", "hubbard_su", "j1j2_su"], Ref("index.md"))
-examples_partition_functions =
-    joinpath.(
-        ["2d_ising_partition_function", "3d_ising_partition_function"], Ref("index.md")
-    )
+examples_optimization = joinpath.(
+    ["heisenberg", "bose_hubbard", "xxz", "fermi_hubbard"], Ref("index.md")
+)
+examples_time_evolution = joinpath.(
+    ["heisenberg_su", "hubbard_su", "j1j2_su"], Ref("index.md")
+)
+examples_partition_functions = joinpath.(
+    ["2d_ising_partition_function", "3d_ising_partition_function"], Ref("index.md")
+)
 examples_boundary_mps = joinpath.(["boundary_mps"], Ref("index.md"))
 
 makedocs(;

src/PEPSKit.jl

@@ -4,6 +4,7 @@ using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf
 using Compat
 using Accessors: @set, @reset
 using VectorInterface
+import VectorInterface as VI
 using TensorKit, KrylovKit, OptimKit, TensorOperations
 using ChainRulesCore, Zygote
 using LoggingExtras

src/algorithms/time_evolution/evoltools.jl

@@ -1,11 +1,11 @@
 """
-    get_expham(dt::Number, H::LocalOperator)
+    get_expham(H::LocalOperator, dt::Number)
 
 Compute `exp(-dt * op)` for each term `op` in `H`,
 and combine them into a new LocalOperator.
 Each `op` in `H` must be a single `TensorMap`.
 """
-function get_expham(dt::Number, H::LocalOperator)
+function get_expham(H::LocalOperator, dt::Number)
     return LocalOperator(
         physicalspace(H), (sites => exp(-dt * op) for (sites, op) in H.terms)...
     )
@@ -181,27 +181,38 @@ Obtain the 3-site gate MPO on the southeast cluster at position `[row, col]`
         c      c+1
 ```
 """
-function _get_gatempo_se(gate::LocalOperator, row::Int, col::Int)
-    Nr, Nc = size(gate.lattice)
+function _get_gatempo_se(ham::LocalOperator, dt::Number, row::Int, col::Int)
+    Nr, Nc = size(ham)
     @assert 1 <= row <= Nr && 1 <= col <= Nc
-    unit = id(space(gate.terms[1].second, 1))
-    sites = (
+    sites = [
         CartesianIndex(row, col),
         CartesianIndex(row, col + 1),
         CartesianIndex(row - 1, col + 1),
-    )
-    nb1x = get_gateterm(gate, (sites[1], sites[2]))
-    nb1y = get_gateterm(gate, (sites[2], sites[3]))
-    nb2 = get_gateterm(gate, (sites[1], sites[3]))
-    op = (1 / 2) * (nb1x ⊗ unit + unit ⊗ nb1y) + permute(nb2 ⊗ unit, ((1, 3, 2), (4, 6, 5)))
+    ]
+    nb1x = get_gateterm(ham, (sites[1], sites[2]))
+    nb1y = get_gateterm(ham, (sites[2], sites[3]))
+    nb2 = get_gateterm(ham, (sites[1], sites[3]))
+    # identity operator at each site
+    units = map(sites) do site
+        site_ = CartesianIndex(mod1(site[1], Nr), mod1(site[2], Nc))
+        return id(physicalspace(ham)[site_])
+    end
+    # when iterating through ┘, └, ┌, ┐ clusters in the unit cell,
+    # NN / NNN bonds are counted 4 / 2 times, respectively.
+    @tensor Odt[i' j' k'; i j k] :=
+        -dt * (
+            (nb1x[i' j'; i j] * units[3][k' k] + units[1][i'; i] * nb1y[j' k'; j k]) / 4 +
+            (nb2[i' k'; i k] * units[2][j'; j]) / 2
+        )
+    op = exp(Odt)
     return gate_to_mpo3(op)
 end
 
 """
 Construct the 3-site gate MPOs on the southeast cluster 
 for 3-site simple update on square lattice.
 """
-function _get_gatempos_se(gate::LocalOperator)
-    Nr, Nc = size(gate.lattice)
-    return collect(_get_gatempo_se(gate, r, c) for r in 1:Nr, c in 1:Nc)
+function _get_gatempos_se(ham::LocalOperator, dt::Number)
+    Nr, Nc = size(ham.lattice)
+    return collect(_get_gatempo_se(ham, dt, r, c) for r in 1:Nr, c in 1:Nc)
 end

src/algorithms/time_evolution/simpleupdate.jl

@@ -9,7 +9,7 @@ Each SU run is converged when the singular value difference becomes smaller than
 $(TYPEDFIELDS)
 """
 struct SimpleUpdate
-    dt::Float64
+    dt::Number
     tol::Float64
     maxiter::Int
     trscheme::TensorKit.TruncationScheme
@@ -29,7 +29,7 @@ Simple update of the x-bond `peps.weights[1,r,c]`.
                 [2,r+1,c]           [2,r+1,c+1]
 ```
 """
-function _su_bondx!(
+function _su_xbond!(
     row::Int,
     col::Int,
     gate::AbstractTensorMap{T,S,2,2},
@@ -94,7 +94,7 @@ function su_iter(
                     direction == 1 ? gate : gate_mirrored,
                     (CartesianIndex(r, 1), CartesianIndex(r, 2)),
                 )
-                ϵ = _su_bondx!(r, 1, term, peps2, alg)
+                ϵ = _su_xbond!(r, 1, term, peps2, alg)
                 peps2.vertices[rp1, 2] = deepcopy(peps2.vertices[r, 1])
                 peps2.vertices[rp1, 1] = deepcopy(peps2.vertices[r, 2])
                 peps2.weights[1, rp1, 2] = deepcopy(peps2.weights[1, r, 1])
@@ -106,7 +106,7 @@ function su_iter(
                     direction == 1 ? gate : gate_mirrored,
                     (CartesianIndex(r, c), CartesianIndex(r, c + 1)),
                 )
-                ϵ = _su_bondx!(r, c, term, peps2, alg)
+                ϵ = _su_xbond!(r, c, term, peps2, alg)
             end
         end
         if direction == 2
@@ -128,7 +128,7 @@ function _simpleupdate2site(
 )
     time_start = time()
     # exponentiating the 2-site Hamiltonian gate
-    gate = get_expham(alg.dt, ham)
+    gate = get_expham(ham, alg.dt)
     wtdiff = 1.0
     wts0 = deepcopy(peps.weights)
     for count in 1:(alg.maxiter)

src/algorithms/time_evolution/simpleupdate3site.jl

@@ -414,13 +414,13 @@ function su3site_iter(
         ),
     )
     peps2 = deepcopy(peps)
-    for i in 1:2
+    for i in 1:4
         for site in CartesianIndices(peps2.vertices)
             r, c = site[1], site[2]
             gs = gatempos[i][r, c]
             _su3site_se!(r, c, gs, peps2, alg)
         end
-        peps2 = (i == 1) ? rotl90(peps2) : rotr90(peps2)
+        peps2 = rotl90(peps2)
     end
     return peps2
 end
@@ -432,9 +432,13 @@ function _simpleupdate3site(
     peps::InfiniteWeightPEPS, ham::LocalOperator, alg::SimpleUpdate; check_interval::Int=500
 )
     time_start = time()
-    gate = get_expham(alg.dt, ham)
-    # convert gates to 3-site MPOs
-    gatempos = [_get_gatempos_se(gate), _get_gatempos_se(rotl90(gate))]
+    # convert Hamiltonian to 3-site exponentiated gate MPOs
+    gatempos = [
+        _get_gatempos_se(ham, alg.dt),
+        _get_gatempos_se(rotl90(ham), alg.dt),
+        _get_gatempos_se(rot180(ham), alg.dt),
+        _get_gatempos_se(rotr90(ham), alg.dt),
+    ]
     wtdiff = 1.0
     wts0 = deepcopy(peps.weights)
     for count in 1:(alg.maxiter)

src/algorithms/toolbox.jl

@@ -272,8 +272,9 @@ passed through to `MPSKit.transfer_spectrum` (e.g. allowing to target the correl
 in a specific symmetry sector).
 
 """
-MPSKit.correlation_length(state, env::CTMRGEnv; num_vals=2, kwargs...) =
+function MPSKit.correlation_length(state, env::CTMRGEnv; num_vals=2, kwargs...)
     _correlation_length(env; num_vals, kwargs...)
+end
 
 function _correlation_length(
     env::CTMRGEnv; num_vals=2, sector=one(sectortype(env)), kwargs...

src/environments/ctmrg_environments.jl

@@ -488,8 +488,6 @@ end
 # big vector. In other words, the associated vector space is not the natural one associated
 # to the original (physical) system, and addition, scaling, etc. are performed element-wise.
 
-import VectorInterface as VI
-
 function VI.scalartype(::Type{CTMRGEnv{C,T}}) where {C,T}
     S₁ = scalartype(C)
     S₂ = scalartype(T)

src/operators/infinitepepo.jl

@@ -177,12 +177,8 @@ end
 
 ## Vector interface
 
-function VectorInterface.scalartype(::Type{NT}) where {NT<:InfinitePEPO}
-    return scalartype(eltype(NT))
-end
-function VectorInterface.zerovector(A::InfinitePEPO)
-    return InfinitePEPO(zerovector(unitcell(A)))
-end
+VI.scalartype(::Type{NT}) where {NT<:InfinitePEPO} = scalartype(eltype(NT))
+VI.zerovector(A::InfinitePEPO) = InfinitePEPO(zerovector(unitcell(A)))
 
 ## (Approximate) equality
 function Base.:(==)(A₁::InfinitePEPO, A₂::InfinitePEPO)

src/operators/localoperator.jl

@@ -105,6 +105,8 @@ function physicalspace(O::LocalOperator)
     return O.lattice
 end
 
+Base.size(O::LocalOperator) = size(physicalspace(O))
+
 # Real and imaginary part
 # -----------------------
 function Base.real(O::LocalOperator)