[ITensors] Pass eigen_perturbation properly to factorize_eigen

src/mps/dmrg.jl

@@ -37,21 +37,20 @@ function dmrg(Hs::Vector{MPO}, psi0::MPS, sweeps::Sweeps; kwargs...)
   return dmrg(PHS, psi0, sweeps; kwargs...)
 end
 
-function dmrg(H::MPO, Ms::Vector{MPS}, psi0::MPS, sweeps::Sweeps; kwargs...)
+function dmrg(H::MPO, Ms::Vector{MPS}, psi0::MPS, sweeps::Sweeps; weight=1.0, kwargs...)
   check_hascommoninds(siteinds, H, psi0)
   check_hascommoninds(siteinds, H, psi0')
   for M in Ms
     check_hascommoninds(siteinds, M, psi0)
   end
   H = permute(H, (linkind, siteinds, linkind))
   Ms .= permute.(Ms, Ref((linkind, siteinds, linkind)))
-  weight = get(kwargs, :weight, 1.0)
   if weight <= 0.0
     error(
       "weight parameter should be > 0.0 in call to excited-state dmrg (value passed was weight=$weight)",
     )
   end
-  PMM = ProjMPO_MPS(H, Ms; weight=weight)
+  PMM = ProjMPO_MPS(H, Ms; weight)
   return dmrg(PMM, psi0, sweeps; kwargs...)
 end
 
@@ -154,7 +153,24 @@ Optional keyword arguments:
     See the `KrylovKit.jl` documention on the `eigsolve` function for more details:
     [KrylovKit.eigsolve](https://jutho.github.io/KrylovKit.jl/stable/man/eig/#KrylovKit.eigsolve).
 """
-function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)
+function dmrg(
+  PH,
+  psi0::MPS,
+  sweeps::Sweeps;
+  which_decomp=nothing,
+  svd_alg=nothing,
+  obs=NoObserver(),
+  outputlevel=1,
+  write_when_maxdim_exceeds=nothing,
+  write_path=tempdir(),
+  # eigsolve kwargs
+  eigsolve_tol=1e-14,
+  eigsolve_krylovdim=3,
+  eigsolve_maxiter=1,
+  eigsolve_verbosity=0,
+  eigsolve_which_eigenvalue=:SR,
+  ishermitian=true,
+)
   if length(psi0) == 1
     error(
       "`dmrg` currently does not support system sizes of 1. You can diagonalize the MPO tensor directly with tools like `LinearAlgebra.eigen`, `KrylovKit.eigsolve`, etc.",
@@ -168,52 +184,6 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)
     checkflux(PH)
   end
 
-  which_decomp::Union{String,Nothing} = get(kwargs, :which_decomp, nothing)
-  svd_alg::String = get(kwargs, :svd_alg, "divide_and_conquer")
-  obs = get(kwargs, :observer, NoObserver())
-  outputlevel::Int = get(kwargs, :outputlevel, 1)
-
-  write_when_maxdim_exceeds::Union{Int,Nothing} = get(
-    kwargs, :write_when_maxdim_exceeds, nothing
-  )
-  write_path = get(kwargs, :write_path, tempdir())
-
-  # eigsolve kwargs
-  eigsolve_tol::Number = get(kwargs, :eigsolve_tol, 1e-14)
-  eigsolve_krylovdim::Int = get(kwargs, :eigsolve_krylovdim, 3)
-  eigsolve_maxiter::Int = get(kwargs, :eigsolve_maxiter, 1)
-  eigsolve_verbosity::Int = get(kwargs, :eigsolve_verbosity, 0)
-
-  ishermitian::Bool = get(kwargs, :ishermitian, true)
-
-  # TODO: add support for targeting other states with DMRG
-  # (such as the state with the largest eigenvalue)
-  # get(kwargs, :eigsolve_which_eigenvalue, :SR)
-  eigsolve_which_eigenvalue::Symbol = :SR
-
-  # TODO: use this as preferred syntax for passing arguments
-  # to eigsolve
-  #default_eigsolve_args = (tol = 1e-14, krylovdim = 3, maxiter = 1,
-  #                         verbosity = 0, ishermitian = true,
-  #                         which_eigenvalue = :SR)
-  #eigsolve = get(kwargs, :eigsolve, default_eigsolve_args)
-
-  # Keyword argument deprecations
-  if haskey(kwargs, :maxiter)
-    error("""maxiter keyword has been replaced by eigsolve_krylovdim.
-             Note: compared to the C++ version of ITensor,
-             setting eigsolve_krylovdim 3 is the same as setting
-             a maxiter of 2.""")
-  end
-
-  if haskey(kwargs, :errgoal)
-    error("errgoal keyword has been replaced by eigsolve_tol.")
-  end
-
-  if haskey(kwargs, :quiet)
-    error("quiet keyword has been replaced by outputlevel")
-  end
-
   psi = copy(psi0)
   N = length(psi)
   if !isortho(psi) || orthocenter(psi) != 1
@@ -341,15 +311,15 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)
         sweep_is_done = (b == 1 && ha == 2)
         measure!(
           obs;
-          energy=energy,
-          psi=psi,
+          energy,
+          psi,
           projected_operator=PH,
           bond=b,
           sweep=sw,
           half_sweep=ha,
-          spec=spec,
-          outputlevel=outputlevel,
-          sweep_is_done=sweep_is_done,
+          spec,
+          outputlevel,
+          sweep_is_done,
         )
       end
     end

src/tensor_operations/matrix_decomposition.jl

@@ -646,7 +646,6 @@ function factorize_eigen(
   cutoff=nothing,
   tags=nothing,
 )
-  delta_A2 = eigen_perturbation
   if ortho == "left"
     Lis = commoninds(A, indices(Linds...))
   elseif ortho == "right"
@@ -657,11 +656,11 @@ function factorize_eigen(
   end
   simLis = sim(Lis)
   A2 = A * replaceinds(dag(A), Lis, simLis)
-  if !isnothing(delta_A2)
+  if !isnothing(eigen_perturbation)
     # This assumes delta_A2 has indices:
     # (Lis..., prime(Lis)...)
-    delta_A2 = replaceinds(delta_A2, Lis, dag(simLis))
-    noprime!(delta_A2)
+    delta_A2 = replaceinds(eigen_perturbation, Lis, dag(simLis))
+    delta_A2 = noprime(delta_A2)
     A2 += delta_A2
   end
   F = eigen(A2, Lis, simLis; ishermitian=true, mindim, maxdim, cutoff, tags)
@@ -782,7 +781,9 @@ function factorize(
     end
     L, R, spec = LR
   elseif which_decomp == "eigen"
-    L, R, spec = factorize_eigen(A, Linds...; mindim, maxdim, cutoff, tags, ortho)
+    L, R, spec = factorize_eigen(
+      A, Linds...; mindim, maxdim, cutoff, tags, ortho, eigen_perturbation
+    )
   elseif which_decomp == "qr"
     L, R = factorize_qr(A, Linds...; ortho, tags)
     spec = Spectrum(nothing, 0.0)