From 0fc3cfbb1d1ef15042dea8e809384185c34fc2ea Mon Sep 17 00:00:00 2001 From: mtfishman Date: Mon, 6 Nov 2023 17:51:44 -0500 Subject: [PATCH] [ITensors] Pass eigen_perturbation properly to factorize_eigen --- src/mps/dmrg.jl | 80 ++++++------------- src/tensor_operations/matrix_decomposition.jl | 11 +-- 2 files changed, 31 insertions(+), 60 deletions(-) diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl index e979b8bce8..d74c23bdef 100644 --- a/src/mps/dmrg.jl +++ b/src/mps/dmrg.jl @@ -37,7 +37,7 @@ 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 @@ -45,13 +45,12 @@ function dmrg(H::MPO, Ms::Vector{MPS}, psi0::MPS, sweeps::Sweeps; kwargs...) 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 diff --git a/src/tensor_operations/matrix_decomposition.jl b/src/tensor_operations/matrix_decomposition.jl index 378bdfa826..6f9b6d67c5 100644 --- a/src/tensor_operations/matrix_decomposition.jl +++ b/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)