[ITensorMPS] Move ITensorTDVP code into ITensorMPS module (#1366)

Project.toml

@@ -17,6 +17,7 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
+Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0"
 PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

src/ITensorMPS/ITensorMPS.jl

@@ -1,6 +1,7 @@
 module ITensorMPS
 
-include("usings.jl")
+using ..ITensors
+
 include("imports.jl")
 include("exports.jl")
 include("abstractmps.jl")
@@ -23,4 +24,19 @@ include("autompo/opsum_to_mpo.jl")
 include("autompo/opsum_to_mpo_qn.jl")
 include("deprecated.jl")
 
+include("abstractprojmpo/projmpo_apply.jl")
+include("abstractprojmpo/projmps2.jl")
+include("abstractprojmpo/projmpo_mps2.jl")
+include("defaults.jl")
+include("update_observer.jl")
+include("solver_utils.jl")
+include("tdvporder.jl")
+include("tdvpinfo.jl")
+include("sweep_update.jl")
+include("alternating_update.jl")
+include("tdvp.jl")
+include("dmrg_x.jl")
+include("contract_mpo_mps.jl")
+include("linsolve.jl")
+
 end # module ITensorMPS

src/ITensorMPS/abstractmps.jl

@@ -1,3 +1,6 @@
+using IsApprox: Approx, IsApprox
+using NDTensors: using_auto_fermion, scalartype, tensor
+
 abstract type AbstractMPS end
 
 """
@@ -45,9 +48,9 @@ have type `ComplexF64`, return `ComplexF64`.
 """
 promote_itensor_eltype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m)
 
-scalartype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m)
-scalartype(m::Array{ITensor}) = LinearAlgebra.promote_leaf_eltypes(m)
-scalartype(m::Array{<:Array{ITensor}}) = LinearAlgebra.promote_leaf_eltypes(m)
+NDTensors.scalartype(m::AbstractMPS) = LinearAlgebra.promote_leaf_eltypes(m)
+NDTensors.scalartype(m::Array{ITensor}) = LinearAlgebra.promote_leaf_eltypes(m)
+NDTensors.scalartype(m::Array{<:Array{ITensor}}) = LinearAlgebra.promote_leaf_eltypes(m)
 
 """
     eltype(m::MPS)
@@ -1141,7 +1144,7 @@ Same as [`inner`](@ref).
 
 See also [`loginner`](@ref), [`logdot`](@ref).
 """
-function dot(M1::MPST, M2::MPST; kwargs...) where {MPST<:AbstractMPS}
+function LinearAlgebra.dot(M1::MPST, M2::MPST; kwargs...) where {MPST<:AbstractMPS}
   return _log_or_not_dot(M1, M2, false; kwargs...)
 end
 
@@ -1319,7 +1322,7 @@ lognorm_ψ[1] == -Inf # There was an infinite norm
 
 See also [`normalize`](@ref), [`norm`](@ref), [`lognorm`](@ref).
 """
-function normalize!(M::AbstractMPS; (lognorm!)=[])
+function LinearAlgebra.normalize!(M::AbstractMPS; (lognorm!)=[])
   c = ortho_lims(M)
   lognorm_M = lognorm(M)
   push!(lognorm!, lognorm_M)

src/ITensorMPS/abstractprojmpo/projmpo_apply.jl

@@ -0,0 +1,89 @@
+using ITensors: ITensor
+
+"""
+A ProjMPOApply represents the application of an
+MPO `H` onto an MPS `psi0` but "projected" by
+the basis of a different MPS `psi` (which
+could be an approximation to H|psi>).
+
+As an implementation of the AbstractProjMPO
+type, it supports multiple `nsite` values for
+one- and two-site algorithms.
+
+```
+     *--*--*-      -*--*--*--*--*--* <psi|
+     |  |  |  |  |  |  |  |  |  |  |
+     h--h--h--h--h--h--h--h--h--h--h H  
+     |  |  |  |  |  |  |  |  |  |  |
+     o--o--o-      -o--o--o--o--o--o |psi0>
+```
+"""
+mutable struct ProjMPOApply <: AbstractProjMPO
+  lpos::Int
+  rpos::Int
+  nsite::Int
+  psi0::MPS
+  H::MPO
+  LR::Vector{ITensor}
+end
+
+function ProjMPOApply(psi0::MPS, H::MPO)
+  return ProjMPOApply(0, length(H) + 1, 2, psi0, H, Vector{ITensor}(undef, length(H)))
+end
+
+function Base.copy(P::ProjMPOApply)
+  return ProjMPOApply(P.lpos, P.rpos, P.nsite, copy(P.psi0), copy(P.H), copy(P.LR))
+end
+
+function set_nsite!(P::ProjMPOApply, nsite)
+  P.nsite = nsite
+  return P
+end
+
+function makeL!(P::ProjMPOApply, psi::MPS, k::Int)
+  # Save the last `L` that is made to help with caching
+  # for DiskProjMPO
+  ll = P.lpos
+  if ll ≥ k
+    # Special case when nothing has to be done.
+    # Still need to change the position if lproj is
+    # being moved backward.
+    P.lpos = k
+    return nothing
+  end
+  # Make sure ll is at least 0 for the generic logic below
+  ll = max(ll, 0)
+  L = lproj(P)
+  while ll < k
+    L = L * P.psi0[ll + 1] * P.H[ll + 1] * dag(psi[ll + 1])
+    P.LR[ll + 1] = L
+    ll += 1
+  end
+  # Needed when moving lproj backward.
+  P.lpos = k
+  return P
+end
+
+function makeR!(P::ProjMPOApply, psi::MPS, k::Int)
+  # Save the last `R` that is made to help with caching
+  # for DiskProjMPO
+  rl = P.rpos
+  if rl ≤ k
+    # Special case when nothing has to be done.
+    # Still need to change the position if rproj is
+    # being moved backward.
+    P.rpos = k
+    return nothing
+  end
+  N = length(P.H)
+  # Make sure rl is no bigger than `N + 1` for the generic logic below
+  rl = min(rl, N + 1)
+  R = rproj(P)
+  while rl > k
+    R = R * P.psi0[rl - 1] * P.H[rl - 1] * dag(psi[rl - 1])
+    P.LR[rl - 1] = R
+    rl -= 1
+  end
+  P.rpos = k
+  return P
+end

src/ITensorMPS/abstractprojmpo/projmpo_mps2.jl

@@ -0,0 +1,46 @@
+using ITensors: ITensors, contract
+
+mutable struct ProjMPO_MPS2 <: AbstractProjMPO
+  PH::ProjMPO
+  Ms::Vector{ProjMPS2}
+end
+
+function ProjMPO_MPS2(H::MPO, M::MPS)
+  return ProjMPO_MPS2(ProjMPO(H), [ProjMPS2(M)])
+end
+
+function ProjMPO_MPS2(H::MPO, Mv::Vector{MPS})
+  return ProjMPO_MPS2(ProjMPO(H), [ProjMPS2(m) for m in Mv])
+end
+
+Base.copy(P::ProjMPO_MPS2) = ProjMPO_MPS2(copy(P.PH), copy(P.Ms))
+
+nsite(P::ProjMPO_MPS2) = nsite(P.PH)
+
+function set_nsite!(P::ProjMPO_MPS2, nsite)
+  set_nsite!(P.PH, nsite)
+  for m in P.Ms
+    set_nsite!(m, nsite)
+  end
+  return P
+end
+
+function makeL!(P::ProjMPO_MPS2, psi::MPS, k::Int)
+  makeL!(P.PH, psi, k)
+  for m in P.Ms
+    makeL!(m, psi, k)
+  end
+  return P
+end
+
+function makeR!(P::ProjMPO_MPS2, psi::MPS, k::Int)
+  makeR!(P.PH, psi, k)
+  for m in P.Ms
+    makeR!(m, psi, k)
+  end
+  return P
+end
+
+ITensors.contract(P::ProjMPO_MPS2, v::ITensor) = contract(P.PH, v)
+
+proj_mps(P::ProjMPO_MPS2) = [proj_mps(m) for m in P.Ms]

src/ITensorMPS/abstractprojmpo/projmps2.jl

@@ -0,0 +1,122 @@
+using ITensors: ITensors, ITensor, dag, dim, prime
+
+"""
+Holds the following data where psi
+is the MPS being optimized and M is the 
+MPS held constant by the ProjMPS.
+```
+     o--o--o--o--o--o--o--o--o--o--o <M|
+     |  |  |  |  |  |  |  |  |  |  |
+     *--*--*-      -*--*--*--*--*--* |psi>
+```
+"""
+mutable struct ProjMPS2 <: AbstractProjMPO
+  lpos::Int
+  rpos::Int
+  nsite::Int
+  M::MPS
+  LR::Vector{ITensor}
+end
+
+function ProjMPS2(M::MPS)
+  return ProjMPS2(0, length(M) + 1, 2, M, Vector{ITensor}(undef, length(M)))
+end
+
+Base.length(P::ProjMPS2) = length(P.M)
+
+function Base.copy(P::ProjMPS2)
+  return ProjMPS2(P.lpos, P.rpos, P.nsite, copy(P.M), copy(P.LR))
+end
+
+function set_nsite!(P::ProjMPS2, nsite)
+  P.nsite = nsite
+  return P
+end
+
+function makeL!(P::ProjMPS2, psi::MPS, k::Int)
+  # Save the last `L` that is made to help with caching
+  # for DiskProjMPO
+  ll = P.lpos
+  if ll ≥ k
+    # Special case when nothing has to be done.
+    # Still need to change the position if lproj is
+    # being moved backward.
+    P.lpos = k
+    return nothing
+  end
+  # Make sure ll is at least 0 for the generic logic below
+  ll = max(ll, 0)
+  L = lproj(P)
+  while ll < k
+    L = L * psi[ll + 1] * dag(prime(P.M[ll + 1], "Link"))
+    P.LR[ll + 1] = L
+    ll += 1
+  end
+  # Needed when moving lproj backward.
+  P.lpos = k
+  return P
+end
+
+function makeR!(P::ProjMPS2, psi::MPS, k::Int)
+  # Save the last `R` that is made to help with caching
+  # for DiskProjMPO
+  rl = P.rpos
+  if rl ≤ k
+    # Special case when nothing has to be done.
+    # Still need to change the position if rproj is
+    # being moved backward.
+    P.rpos = k
+    return nothing
+  end
+  N = length(P.M)
+  # Make sure rl is no bigger than `N + 1` for the generic logic below
+  rl = min(rl, N + 1)
+  R = rproj(P)
+  while rl > k
+    R = R * psi[rl - 1] * dag(prime(P.M[rl - 1], "Link"))
+    P.LR[rl - 1] = R
+    rl -= 1
+  end
+  P.rpos = k
+  return P
+end
+
+function ITensors.contract(P::ProjMPS2, v::ITensor)
+  itensor_map = Union{ITensor,OneITensor}[lproj(P)]
+  append!(itensor_map, [prime(t, "Link") for t in P.M[site_range(P)]])
+  push!(itensor_map, rproj(P))
+
+  # Reverse the contraction order of the map if
+  # the first tensor is a scalar (for example we
+  # are at the left edge of the system)
+  if dim(first(itensor_map)) == 1
+    reverse!(itensor_map)
+  end
+
+  # Apply the map
+  Mv = v
+  for it in itensor_map
+    Mv *= it
+  end
+  return Mv
+end
+
+function proj_mps(P::ProjMPS2)
+  itensor_map = Union{ITensor,OneITensor}[lproj(P)]
+  append!(itensor_map, [dag(prime(t, "Link")) for t in P.M[site_range(P)]])
+  push!(itensor_map, rproj(P))
+
+  # Reverse the contraction order of the map if
+  # the first tensor is a scalar (for example we
+  # are at the left edge of the system)
+  if dim(first(itensor_map)) == 1
+    reverse!(itensor_map)
+  end
+
+  # Apply the map
+  m = ITensor(true)
+  for it in itensor_map
+    m *= it
+  end
+  return m
+end

src/ITensorMPS/adapt.jl

@@ -1 +1,2 @@
-adapt_structure(to, x::Union{MPS,MPO}) = map(xᵢ -> adapt(to, xᵢ), x)
+using Adapt: Adapt
+Adapt.adapt_structure(to, x::Union{MPS,MPO}) = map(xᵢ -> adapt(to, xᵢ), x)