diff --git a/Project.toml b/Project.toml index 7a08250..3064cac 100644 --- a/Project.toml +++ b/Project.toml @@ -6,7 +6,6 @@ version = "0.3.0" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" -ITensorTDVP = "25707e16-a4db-4a07-99d9-4d67b7af0342" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" IsApprox = "28f27b66-4bd8-47e7-9110-e2746eb8bed7" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" @@ -18,10 +17,15 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SerializedElementArrays = "d3ce8812-9567-47e9-a7b5-65a6d70a3065" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" +[weakdeps] +ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" + +[extensions] +ITensorMPSChainRulesCoreExt = "ChainRulesCore" + [compat] Adapt = "4.1.0" Compat = "4.16.0" -ITensorTDVP = "0.4.1" ITensors = "0.7" IsApprox = "2.0.0" KrylovKit = "0.8.1" diff --git a/examples/solvers/Project.toml b/examples/solvers/Project.toml index 4d3b745..ae4dddd 100644 --- a/examples/solvers/Project.toml +++ b/examples/solvers/Project.toml @@ -6,6 +6,3 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0" OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" - -[compat] -ITensors = "0.6.7" diff --git a/ext/ITensorMPSChainRulesCoreExt/abstractmps.jl b/ext/ITensorMPSChainRulesCoreExt/abstractmps.jl index ee6295a..53c6980 100644 --- a/ext/ITensorMPSChainRulesCoreExt/abstractmps.jl +++ b/ext/ITensorMPSChainRulesCoreExt/abstractmps.jl @@ -2,7 +2,7 @@ using Adapt: adapt using ChainRulesCore: ChainRulesCore, HasReverseMode, NoTangent, RuleConfig, rrule_via_ad using ITensors: ITensors, ITensor, dag, hassameinds, inds, itensor, mapprime, replaceprime, swapprime -using ITensors.ITensorMPS: ITensorMPS, MPO, MPS, apply, inner, siteinds +using ITensorMPS: ITensorMPS, MPO, MPS, apply, inner, siteinds using NDTensors: datatype function ChainRulesCore.rrule( diff --git a/ext/ITensorMPSChainRulesCoreExt/indexset.jl b/ext/ITensorMPSChainRulesCoreExt/indexset.jl index a5f4409..154f3f3 100644 --- a/ext/ITensorMPSChainRulesCoreExt/indexset.jl +++ b/ext/ITensorMPSChainRulesCoreExt/indexset.jl @@ -10,7 +10,7 @@ using ITensors: replacetags, setprime, settags -using ITensors.ITensorMPS: MPO, MPS +using ITensorMPS: MPO, MPS for fname in ( :prime, :setprime, :noprime, :replaceprime, :addtags, :removetags, :replacetags, :settags diff --git a/ext/ITensorMPSChainRulesCoreExt/mpo.jl b/ext/ITensorMPSChainRulesCoreExt/mpo.jl index 14cdd6c..3da69ac 100644 --- a/ext/ITensorMPSChainRulesCoreExt/mpo.jl +++ b/ext/ITensorMPSChainRulesCoreExt/mpo.jl @@ -1,6 +1,6 @@ using ChainRulesCore: ChainRulesCore, NoTangent using ITensors: Algorithm, contract, hassameinds, inner, mapprime -using ITensors.ITensorMPS: MPO, MPS, firstsiteinds, siteinds +using ITensorMPS: MPO, MPS, firstsiteinds, siteinds using LinearAlgebra: tr function ChainRulesCore.rrule( diff --git a/ext/ITensorMPSChainRulesCoreExt/mps.jl b/ext/ITensorMPSChainRulesCoreExt/mps.jl index 43f27a4..e0647d4 100644 --- a/ext/ITensorMPSChainRulesCoreExt/mps.jl +++ b/ext/ITensorMPSChainRulesCoreExt/mps.jl @@ -1,4 +1,4 @@ using ChainRulesCore: @non_differentiable using ITensors: Index -using ITensors.ITensorMPS: MPS +using ITensorMPS: MPS @non_differentiable MPS(::Type{<:Number}, sites::Vector{<:Index}, states_) diff --git a/ext/ITensorMPSObserversExt/ITensorMPSObserversExt.jl b/ext/ITensorMPSObserversExt/ITensorMPSObserversExt.jl new file mode 100644 index 0000000..f353b99 --- /dev/null +++ b/ext/ITensorMPSObserversExt/ITensorMPSObserversExt.jl @@ -0,0 +1,9 @@ +module ITensorMPSObserversExt +using Observers: Observers +using Observers.DataFrames: AbstractDataFrame +using ITensorMPS: ITensorMPS + +function ITensorMPS.update_observer!(observer::AbstractDataFrame; kwargs...) + return Observers.update!(observer; kwargs...) +end +end diff --git a/ext/ITensorMPSPackageCompilerExt/ITensorMPSPackageCompilerExt.jl b/ext/ITensorMPSPackageCompilerExt/ITensorMPSPackageCompilerExt.jl new file mode 100644 index 0000000..3bb7b57 --- /dev/null +++ b/ext/ITensorMPSPackageCompilerExt/ITensorMPSPackageCompilerExt.jl @@ -0,0 +1,3 @@ +module ITensorsPackageCompilerExt +include("compile.jl") +end diff --git a/ext/ITensorMPSPackageCompilerExt/compile.jl b/ext/ITensorMPSPackageCompilerExt/compile.jl new file mode 100644 index 0000000..c14d0e9 --- /dev/null +++ b/ext/ITensorMPSPackageCompilerExt/compile.jl @@ -0,0 +1,26 @@ +using NDTensors: @Algorithm_str +using ITensors: ITensors +using PackageCompiler: PackageCompiler + +function ITensors.compile( + ::Algorithm"PackageCompiler"; + dir::AbstractString=ITensors.default_compile_dir(), + filename::AbstractString=ITensors.default_compile_filename(), +) + if !isdir(dir) + println("""The directory "$dir" doesn't exist yet, creating it now.""") + println() + mkdir(dir) + end + path = joinpath(dir, filename) + println( + """Creating the system image "$path" containing the compiled version of ITensorMPS. This may take a few minutes.""", + ) + PackageCompiler.create_sysimage( + :ITensorMPS; + sysimage_path=path, + precompile_execution_file=joinpath(@__DIR__, "precompile_itensormps.jl"), + ) + println(ITensors.compile_note(; dir, filename)) + return path +end diff --git a/ext/ITensorMPSPackageCompilerExt/precompile_itensormps.jl b/ext/ITensorMPSPackageCompilerExt/precompile_itensormps.jl new file mode 100644 index 0000000..ebcadde --- /dev/null +++ b/ext/ITensorMPSPackageCompilerExt/precompile_itensormps.jl @@ -0,0 +1,28 @@ +using ITensorMPS: MPO, OpSum, dmrg, random_mps, siteinds + +# TODO: This uses all of the tests to make +# precompile statements, but takes a long time +# (e.g. 700 seconds). +# Try again with later versions of PackageCompiler +# +# include(joinpath(joinpath(dirname(dirname(@__DIR__)), +# test"), +# "runtests.jl")) + +function main(; N, dmrg_kwargs) + opsum = OpSum() + for j in 1:(N - 1) + opsum += 0.5, "S+", j, "S-", j + 1 + opsum += 0.5, "S-", j, "S+", j + 1 + opsum += "Sz", j, "Sz", j + 1 + end + for conserve_qns in (false, true) + sites = siteinds("S=1", N; conserve_qns) + H = MPO(opsum, sites) + ψ0 = random_mps(sites, j -> isodd(j) ? "↑" : "↓"; linkdims=2) + dmrg(H, ψ0; outputlevel=0, dmrg_kwargs...) + end + return nothing +end + +main(; N=6, dmrg_kwargs=(; nsweeps=3, maxdim=10, cutoff=1e-13)) diff --git a/src/ITensorMPS.jl b/src/ITensorMPS.jl index 82779d4..30d3778 100644 --- a/src/ITensorMPS.jl +++ b/src/ITensorMPS.jl @@ -25,4 +25,24 @@ include("defaults.jl") include("update_observer.jl") include("lattices/lattices.jl") export Lattice, LatticeBond, square_lattice, triangular_lattice +export TimeDependentSum, dmrg_x, expand, linsolve, tdvp, to_vec +include("solvers/ITensorsExtensions.jl") +using .ITensorsExtensions: to_vec +include("solvers/applyexp.jl") +include("solvers/defaults.jl") +include("solvers/update_observer.jl") +include("solvers/timedependentsum.jl") +include("solvers/tdvporder.jl") +include("solvers/sweep_update.jl") +include("solvers/alternating_update.jl") +include("solvers/tdvp.jl") +include("solvers/dmrg.jl") +include("solvers/dmrg_x.jl") +include("solvers/reducedcontractproblem.jl") +include("solvers/contract.jl") +include("solvers/reducedconstantterm.jl") +include("solvers/reducedlinearproblem.jl") +include("solvers/linsolve.jl") +include("solvers/expand.jl") +include("lib/Experimental/src/Experimental.jl") end diff --git a/src/abstractmps.jl b/src/abstractmps.jl index 3b663af..b3db530 100644 --- a/src/abstractmps.jl +++ b/src/abstractmps.jl @@ -1833,7 +1833,7 @@ function setindex!( # into MPS tensors firstsite = first(r) lastsite = last(r) - @assert firstsite ≤ ITensors.orthocenter(ψ) ≤ lastsite + @assert firstsite ≤ ITensorMPS.orthocenter(ψ) ≤ lastsite @assert firstsite ≤ leftlim(ψ) + 1 @assert rightlim(ψ) - 1 ≤ lastsite diff --git a/src/Experimental.jl b/src/lib/Experimental/src/Experimental.jl similarity index 50% rename from src/Experimental.jl rename to src/lib/Experimental/src/Experimental.jl index 42be6dc..4e6fa93 100644 --- a/src/Experimental.jl +++ b/src/lib/Experimental/src/Experimental.jl @@ -1,3 +1,3 @@ module Experimental -using ITensorTDVP: dmrg +include("dmrg.jl") end diff --git a/src/lib/Experimental/src/dmrg.jl b/src/lib/Experimental/src/dmrg.jl new file mode 100644 index 0000000..8796936 --- /dev/null +++ b/src/lib/Experimental/src/dmrg.jl @@ -0,0 +1,17 @@ +using ..ITensorMPS: + MPS, + alternating_update, + compose_observers, + default_observer, + eigsolve_updater, + values_observer + +function dmrg( + operator, init::MPS; updater=eigsolve_updater, (observer!)=default_observer(), kwargs... +) + info_ref! = Ref{Any}() + info_observer! = values_observer(; info=info_ref!) + observer! = compose_observers(observer!, info_observer!) + state = alternating_update(operator, init; updater, observer!, kwargs...) + return info_ref![].eigval, state +end diff --git a/src/solvers/ITensorsExtensions.jl b/src/solvers/ITensorsExtensions.jl new file mode 100644 index 0000000..90ed8b1 --- /dev/null +++ b/src/solvers/ITensorsExtensions.jl @@ -0,0 +1,9 @@ +module ITensorsExtensions +using ITensors: ITensor, array, inds, itensor +function to_vec(x::ITensor) + function to_itensor(x_vec) + return itensor(x_vec, inds(x)) + end + return vec(array(x)), to_itensor +end +end diff --git a/src/solvers/alternating_update.jl b/src/solvers/alternating_update.jl new file mode 100644 index 0000000..75ecc4b --- /dev/null +++ b/src/solvers/alternating_update.jl @@ -0,0 +1,117 @@ +using ITensors: ITensors, permute + +function _extend_sweeps_param(param, nsweeps) + if param isa Number + eparam = fill(param, nsweeps) + else + length(param) == nsweeps && return param + eparam = Vector(undef, nsweeps) + eparam[1:length(param)] = param + eparam[(length(param) + 1):end] .= param[end] + end + return eparam +end + +function process_sweeps(; nsweeps, maxdim, mindim, cutoff, noise) + maxdim = _extend_sweeps_param(maxdim, nsweeps) + mindim = _extend_sweeps_param(mindim, nsweeps) + cutoff = _extend_sweeps_param(cutoff, nsweeps) + noise = _extend_sweeps_param(noise, nsweeps) + return (; maxdim, mindim, cutoff, noise) +end + +function alternating_update( + operator, + init::MPS; + updater, + updater_kwargs=(;), + nsweeps=default_nsweeps(), + checkdone=default_checkdone(), + write_when_maxdim_exceeds=default_write_when_maxdim_exceeds(), + nsite=default_nsite(), + reverse_step=default_reverse_step(), + time_start=default_time_start(), + time_step=default_time_step(), + order=default_order(), + (observer!)=default_observer(), + (sweep_observer!)=default_sweep_observer(), + outputlevel=default_outputlevel(), + normalize=default_normalize(), + maxdim=default_maxdim(), + mindim=default_mindim(), + cutoff=default_cutoff(ITensors.scalartype(init)), + noise=default_noise(), +) + reduced_operator = ITensorMPS.reduced_operator(operator) + if isnothing(nsweeps) + return error("Must specify `nsweeps`.") + end + maxdim, mindim, cutoff, noise = process_sweeps(; nsweeps, maxdim, mindim, cutoff, noise) + forward_order = TDVPOrder(order, Base.Forward) + state = copy(init) + # Keep track of the start of the current time step. + # Helpful for tracking the total time, for example + # when using time-dependent updaters. + # This will be passed as a keyword argument to the + # `updater`. + current_time = time_start + info = nothing + for sweep in 1:nsweeps + if !isnothing(write_when_maxdim_exceeds) && maxdim[sweep] > write_when_maxdim_exceeds + if outputlevel >= 2 + println( + "write_when_maxdim_exceeds = $write_when_maxdim_exceeds and maxdim(sweeps, sw) = $(maxdim(sweeps, sweep)), writing environment tensors to disk", + ) + end + reduced_operator = disk(reduced_operator) + end + sweep_elapsed_time = @elapsed begin + state, reduced_operator, info = sweep_update( + forward_order, + reduced_operator, + state; + updater, + updater_kwargs, + nsite, + current_time, + time_step, + reverse_step, + sweep, + observer!, + normalize, + outputlevel, + maxdim=maxdim[sweep], + mindim=mindim[sweep], + cutoff=cutoff[sweep], + noise=noise[sweep], + ) + end + if !isnothing(time_step) + current_time += time_step + end + update_observer!( + sweep_observer!; state, reduced_operator, sweep, outputlevel, current_time + ) + if outputlevel >= 1 + print("After sweep ", sweep, ":") + print(" maxlinkdim=", maxlinkdim(state)) + @printf(" maxerr=%.2E", info.maxtruncerr) + if !isnothing(current_time) + print(" current_time=", round(current_time; digits=3)) + end + print(" time=", round(sweep_elapsed_time; digits=3)) + println() + flush(stdout) + end + isdone = checkdone(; + state, sweep, outputlevel, observer=observer!, sweep_observer=sweep_observer! + ) + isdone && break + end + return state +end + +# Assume it is already in a reduced basis. +reduced_operator(operator) = operator +reduced_operator(operators::Vector{MPO}) = ProjMPOSum(operators) +reduced_operator(operator::MPO) = ProjMPO(operator) diff --git a/src/solvers/applyexp.jl b/src/solvers/applyexp.jl new file mode 100644 index 0000000..92121b8 --- /dev/null +++ b/src/solvers/applyexp.jl @@ -0,0 +1,131 @@ +using LinearAlgebra: dot +using Printf: @printf + +# +# To Do: +# - implement assembleLanczosVectors +# - check slice ranges - change end value by 1? +# + +function assemble_lanczos_vecs(lanczos_vectors, linear_comb, norm) + #if length(lanczos_vectors) != length(linear_comb) + # @show length(lanczos_vectors) + # @show length(linear_comb) + #end + xt = norm * linear_comb[1] * lanczos_vectors[1] + for i in 2:length(lanczos_vectors) + xt += norm * linear_comb[i] * lanczos_vectors[i] + end + return xt +end + +struct ApplyExpInfo + numops::Int + converged::Int +end + +function applyexp(H, tau::Number, x0; maxiter=30, tol=1e-12, outputlevel=0, normcutoff=1e-7) + # Initialize Lanczos vectors + v1 = copy(x0) + nrm = norm(v1) + v1 /= nrm + lanczos_vectors = [v1] + + ElT = promote_type(typeof(tau), eltype(x0)) + + bigTmat = zeros(ElT, maxiter + 3, maxiter + 3) + + nmatvec = 0 + + v0 = nothing + beta = 0.0 + for iter in 1:maxiter + tmat_size = iter + 1 + + # Matrix-vector multiplication + w = H(v1) + nmatvec += 1 + + avnorm = norm(w) + alpha = dot(w, v1) + + bigTmat[iter, iter] = alpha + + w -= alpha * v1 + if iter > 1 + w -= beta * v0 + end + v0 = copy(v1) + + beta = norm(w) + + # check for Lanczos sequence exhaustion + if abs(beta) < normcutoff + # Assemble the time evolved state + tmat = bigTmat[1:tmat_size, 1:tmat_size] + tmat_exp = exp(tau * tmat) + linear_comb = tmat_exp[:, 1] + xt = assemble_lanczos_vecs(lanczos_vectors, linear_comb, nrm) + return xt, ApplyExpInfo(nmatvec, 1) + end + + # update next lanczos vector + v1 = copy(w) + v1 /= beta + push!(lanczos_vectors, v1) + bigTmat[iter + 1, iter] = beta + bigTmat[iter, iter + 1] = beta + + # Convergence check + if iter > 0 + # Prepare extended T-matrix for exponentiation + tmat_ext_size = tmat_size + 2 + tmat_ext = bigTmat[1:tmat_ext_size, 1:tmat_ext_size] + + tmat_ext[tmat_size - 1, tmat_size] = 0.0 + tmat_ext[tmat_size + 1, tmat_size] = 1.0 + + # Exponentiate extended T-matrix + tmat_ext_exp = exp(tau * tmat_ext) + + ϕ1 = abs(nrm * tmat_ext_exp[tmat_size, 1]) + ϕ2 = abs(nrm * tmat_ext_exp[tmat_size + 1, 1] * avnorm) + + if ϕ1 > 10 * ϕ2 + error = ϕ2 + elseif (ϕ1 > ϕ2) + error = (ϕ1 * ϕ2) / (ϕ1 - ϕ2) + else + error = ϕ1 + end + + if outputlevel >= 3 + @printf(" Iteration: %d, Error: %.2E\n", iter, error) + end + + if ((error < tol) || (iter == maxiter)) + converged = 1 + if (iter == maxiter) + println("warning: applyexp not converged in $maxiter steps") + converged = 0 + end + + # Assemble the time evolved state + linear_comb = tmat_ext_exp[:, 1] + xt = assemble_lanczos_vecs(lanczos_vectors, linear_comb, nrm) + + if outputlevel >= 3 + println(" Number of iterations: $iter") + end + + return xt, ApplyExpInfo(nmatvec, converged) + end + end # end convergence test + end # iter + + if outputlevel >= 0 + println("In applyexp, number of matrix-vector multiplies: ", nmatvec) + end + + return x0 +end diff --git a/src/solvers/contract.jl b/src/solvers/contract.jl new file mode 100644 index 0000000..1919fdd --- /dev/null +++ b/src/solvers/contract.jl @@ -0,0 +1,41 @@ +using ITensors: ITensors, Index, ITensor, @Algorithm_str, commoninds, contract, hasind, sim + +function contract_operator_state_updater(operator, init; internal_kwargs) + # TODO: Use `contract(operator)`. + state = ITensor(true) + for j in (operator.lpos + 1):(operator.rpos - 1) + state *= operator.input_state[j] + end + state = contract(operator, state) + return state, (;) +end + +function default_contract_init(operator::MPO, input_state::MPS) + input_state = deepcopy(input_state) + s = only.(siteinds(uniqueinds, operator, input_state)) + # TODO: Fix issue with `replace_siteinds`, seems to be modifying in-place. + return replace_siteinds(deepcopy(input_state), s) +end + +function ITensors.contract( + ::Algorithm"fit", + operator::MPO, + input_state::MPS; + init=default_contract_init(operator, input_state), + kwargs..., +) + # Fix siteinds of `init` if needed. + # This is needed to work around an issue that `apply` + # can't be customized right now, and just uses the same `init` + # as that of `contract`. + # TODO: Allow customization of `apply` and remove this. + s = only.(siteinds(uniqueinds, operator, input_state)) + if !all(p -> p[1] == [2], zip(s, siteinds(init))) + # TODO: Fix issue with `replace_siteinds`, seems to be modifying in-place. + init = replace_siteinds(deepcopy(init), s) + end + reduced_operator = ReducedContractProblem(input_state, operator) + return alternating_update( + reduced_operator, init; updater=contract_operator_state_updater, kwargs... + ) +end diff --git a/src/solvers/defaults.jl b/src/solvers/defaults.jl new file mode 100644 index 0000000..6211125 --- /dev/null +++ b/src/solvers/defaults.jl @@ -0,0 +1,15 @@ +using Compat: Returns + +default_nsweeps() = nothing +default_checkdone() = Returns(false) +default_write_when_maxdim_exceeds() = nothing +default_nsite() = 2 +default_reverse_step() = false +default_time_start() = nothing +default_time_step() = nothing +default_order() = 2 +default_observer() = EmptyObserver() +default_sweep_observer() = EmptyObserver() +default_outputlevel() = 0 +default_normalize() = false +default_sweep() = 1 diff --git a/src/solvers/dmrg.jl b/src/solvers/dmrg.jl new file mode 100644 index 0000000..90fc9df --- /dev/null +++ b/src/solvers/dmrg.jl @@ -0,0 +1,25 @@ +using ITensors: ITensors +using KrylovKit: eigsolve + +function eigsolve_updater( + operator, + init; + internal_kwargs, + which_eigval=:SR, + ishermitian=true, + tol=10^2 * eps(real(ITensors.scalartype(init))), + krylovdim=3, + maxiter=1, + verbosity=0, + eager=false, +) + howmany = 1 + eigvals, eigvecs, info = eigsolve( + operator, init, howmany, which_eigval; ishermitian, tol, krylovdim, maxiter, verbosity + ) + return eigvecs[1], (; info, eigval=eigvals[1]) +end + +# A version of `dmrg` based on `alternating_update` and `eigsolve_updater` +# is defined in `src/lib/Experimental` as `ITensorMPS.Experimental.dmrg` to not conflict +# with `ITensorMPS.dmrg`. diff --git a/src/solvers/dmrg_x.jl b/src/solvers/dmrg_x.jl new file mode 100644 index 0000000..48201e6 --- /dev/null +++ b/src/solvers/dmrg_x.jl @@ -0,0 +1,23 @@ +using ITensors: array, contract, dag, uniqueind, onehot +using LinearAlgebra: eigen + +function eigen_updater(operator, state; internal_kwargs) + contracted_operator = contract(operator, ITensor(true)) + d, u = eigen(contracted_operator; ishermitian=true) + u_ind = uniqueind(u, contracted_operator) + u′_ind = uniqueind(d, u) + max_overlap, max_index = findmax(abs, array(state * dag(u))) + u_max = u * dag(onehot(eltype(u), u_ind => max_index)) + d_max = d[u′_ind => max_index, u_ind => max_index] + return u_max, (; eigval=d_max) +end + +function dmrg_x( + operator, state::MPS; updater=eigen_updater, (observer!)=default_observer(), kwargs... +) + info_ref = Ref{Any}() + info_observer = values_observer(; info=info_ref) + observer = compose_observers(observer!, info_observer) + eigvec = alternating_update(operator, state; updater, (observer!)=observer, kwargs...) + return info_ref[].eigval, eigvec +end diff --git a/src/solvers/expand.jl b/src/solvers/expand.jl new file mode 100644 index 0000000..c970e45 --- /dev/null +++ b/src/solvers/expand.jl @@ -0,0 +1,161 @@ +using Adapt: adapt +using ITensors: + ITensors, + Algorithm, + Index, + ITensor, + @Algorithm_str, + δ, + commonind, + dag, + denseblocks, + directsum, + hasqns, + prime, + scalartype, + uniqueinds +using LinearAlgebra: normalize, svd, tr +using NDTensors: unwrap_array_type + +# Possible improvements: +# - Allow a maxdim argument to be passed to `expand`. +# - Current behavior is letting bond dimension get too +# big when used in imaginary time evolution. +# - Consider switching the default to variational/fit apply +# when building Krylov vectors. +# - Use (1-tau*operator)|state> to generate Krylov vectors +# instead of operator|state>. Is that needed? + +function expand(state, reference; alg, kwargs...) + return expand(Algorithm(alg), state, reference; kwargs...) +end + +function expand_cutoff_doctring() + return """ + The cutoff is used to control the truncation of the expanded + basis and defaults to half the precision of the scalar type + of the input state, i.e. ~1e-8 for `Float64`. + """ +end + +function expand_warning_doctring() + return """ + !!! warning + Users are not given many customization options just yet as we + gain more experience on the right balance between efficacy of the + expansion and performance in various scenarios, and default values + and keyword arguments are subject to change as we learn more about + how to best use the method. + """ +end + +function expand_citation_docstring() + return """ + [^global_expansion]: Time Dependent Variational Principle with Ancillary Krylov Subspace. Mingru Yang, Steven R. White, [arXiv:2005.06104](https://arxiv.org/abs/2005.06104) + """ +end + +""" + expand(state::MPS, references::Vector{MPS}; alg="orthogonalize", cutoff) + +Given an MPS `state` and a collection of MPS `references`, +returns an MPS which is equal to `state` +(has fidelity 1 with `state`) but whose MPS basis +is expanded to contain a portion of the basis of +the `references` that is orthogonal to the MPS basis of `state`. +See [^global_expansion] for more details. + +$(expand_cutoff_doctring()) + +$(expand_warning_doctring()) + +$(expand_citation_docstring()) +""" +function expand( + ::Algorithm"orthogonalize", + state::MPS, + references::Vector{MPS}; + cutoff=(√(eps(real(scalartype(state))))), +) + n = length(state) + state = orthogonalize(state, n) + references = map(reference -> orthogonalize(reference, n), references) + s = siteinds(state) + for j in reverse(2:n) + # SVD state[j] to compute basisⱼ + linds = [s[j - 1]; linkinds(state, j - 1)] + _, λⱼ, basisⱼ = svd(state[j], linds; righttags="bψ_$j,Link") + rinds = uniqueinds(basisⱼ, λⱼ) + # Make projectorⱼ + idⱼ = prod(rinds) do r + return adapt(unwrap_array_type(basisⱼ), denseblocks(δ(scalartype(state), r', dag(r)))) + end + projectorⱼ = idⱼ - prime(basisⱼ, rinds) * dag(basisⱼ) + # Sum reference density matrices + ρⱼ = sum(reference -> prime(reference[j], rinds) * dag(reference[j]), references) + ρⱼ /= tr(ρⱼ) + # Apply projectorⱼ + ρⱼ_projected = apply(apply(projectorⱼ, ρⱼ), projectorⱼ) + expanded_basisⱼ = basisⱼ + if norm(ρⱼ_projected) > 10^3 * eps(real(scalartype(state))) + # Diagonalize projected density matrix ρⱼ_projected + # to compute reference_basisⱼ, which spans part of right basis + # of references which is orthogonal to right basis of state + dⱼ, reference_basisⱼ = eigen( + ρⱼ_projected; cutoff, ishermitian=true, righttags="bϕ_$j,Link" + ) + state_indⱼ = only(commoninds(basisⱼ, λⱼ)) + reference_indⱼ = only(commoninds(reference_basisⱼ, dⱼ)) + expanded_basisⱼ, expanded_indⱼ = directsum( + basisⱼ => state_indⱼ, reference_basisⱼ => reference_indⱼ + ) + end + # Shift ortho center one site left using dag(expanded_basisⱼ) + # and replace tensor at site j with expanded_basisⱼ + state[j - 1] = state[j - 1] * (state[j] * dag(expanded_basisⱼ)) + state[j] = expanded_basisⱼ + for reference in references + reference[j - 1] = reference[j - 1] * (reference[j] * dag(expanded_basisⱼ)) + reference[j] = expanded_basisⱼ + end + end + return state +end + +""" + expand(state::MPS, reference::MPO; alg="global_krylov", krylovdim=2, cutoff) + +Given an MPS `state` and an MPO `reference`, +returns an MPS which is equal to `state` +(has fidelity 1 with state) but whose MPS basis +is expanded to contain a portion of the basis of +the Krylov space built by repeated applications of +`reference` to `state` that is orthogonal +to the MPS basis of `state`. +The `reference` operator is applied to `state` `krylovdim` +number of times, with a default of 2 which should give +a good balance between reliability and performance. +See [^global_expansion] for more details. + +$(expand_cutoff_doctring()) + +$(expand_warning_doctring()) + +$(expand_citation_docstring()) +""" +function expand( + ::Algorithm"global_krylov", + state::MPS, + operator::MPO; + krylovdim=2, + cutoff=(√(eps(real(scalartype(state))))), + apply_kwargs=(; maxdim=maxlinkdim(state) + 1), +) + # TODO: Try replacing this logic with `Base.accumulate`. + references = Vector{MPS}(undef, krylovdim) + for k in 1:krylovdim + previous_reference = get(references, k - 1, state) + references[k] = normalize(apply(operator, previous_reference; apply_kwargs...)) + end + return expand(state, references; alg="orthogonalize", cutoff) +end diff --git a/src/solvers/linsolve.jl b/src/solvers/linsolve.jl new file mode 100644 index 0000000..9ed68f3 --- /dev/null +++ b/src/solvers/linsolve.jl @@ -0,0 +1,50 @@ +using KrylovKit: KrylovKit, linsolve + +function linsolve_updater(problem, init; internal_kwargs, coefficients, kwargs...) + x, info = linsolve( + operator(problem), + constant_term(problem), + init, + coefficients[1], + coefficients[2]; + kwargs..., + ) + return x, (; info) +end + +""" +Compute a solution x to the linear system: + +(a₀ + a₁ * A)*x = b + +using starting guess x₀. Leaving a₀, a₁ +set to their default values solves the +system A*x = b. + +To adjust the balance between accuracy of solution +and speed of the algorithm, it is recommed to first try +adjusting the updater keyword arguments as descibed below. + +Keyword arguments: + - `nsweeps`, `cutoff`, `maxdim`, etc. (like for other MPO/MPS updaters). + - `updater_kwargs=(;)` - a `NamedTuple` containing keyword arguments that will get forwarded to the local updater, + in this case `KrylovKit.linsolve` which is a GMRES linear updater. For example: + ```juli + linsolve(A, b, x; maxdim=100, cutoff=1e-8, nsweeps=10, updater_kwargs=(; ishermitian=true, tol=1e-6, maxiter=20, krylovdim=30)) + ``` + See `KrylovKit.jl` documentation for more details on available keyword arguments. +""" +function KrylovKit.linsolve( + operator, + constant_term::MPS, + init::MPS, + coefficient1::Number=false, + coefficient2::Number=true; + updater=linsolve_updater, + updater_kwargs=(;), + kwargs..., +) + reduced_problem = ReducedLinearProblem(operator, constant_term) + updater_kwargs = (; coefficients=(coefficient1, coefficient2), updater_kwargs...) + return alternating_update(reduced_problem, init; updater, updater_kwargs, kwargs...) +end diff --git a/src/solvers/reducedconstantterm.jl b/src/solvers/reducedconstantterm.jl new file mode 100644 index 0000000..52c882b --- /dev/null +++ b/src/solvers/reducedconstantterm.jl @@ -0,0 +1,152 @@ +using ITensors: ITensors, ITensor, dag, dim, prime + +""" +Holds the following data where basis +is the MPS being optimized and state is the +MPS held constant by the ProjMPS. +``` + o--o--o--o--o--o--o--o--o--o--o +``` +""" +mutable struct ReducedConstantTerm <: AbstractProjMPO + lpos::Int + rpos::Int + nsite::Int + state::MPS + environments::Vector{ITensor} +end + +function ReducedConstantTerm(state::MPS) + lpos = 0 + rpos = length(state) + 1 + nsite = 2 + environments = Vector{ITensor}(undef, length(state)) + return ReducedConstantTerm(lpos, rpos, nsite, state, environments) +end + +function Base.getproperty(reduced_state::ReducedConstantTerm, sym::Symbol) + # This is for compatibility with `AbstractProjMPO`. + # TODO: Don't use `reduced_state.H`, `reduced_state.LR`, etc. + # in `AbstractProjMPO`. + if sym == :LR + return getfield(reduced_state, :environments) + end + return getfield(reduced_state, sym) +end + +Base.length(reduced_state::ReducedConstantTerm) = length(reduced_state.state) + +function Base.copy(reduced_state::ReducedConstantTerm) + return ReducedConstantTerm( + reduced_state.lpos, + reduced_state.rpos, + reduced_state.nsite, + copy(reduced_state.state), + copy(reduced_state.environments), + ) +end + +function ITensorMPS.set_nsite!(reduced_state::ReducedConstantTerm, nsite) + reduced_state.nsite = nsite + return reduced_state +end + +function ITensorMPS.makeL!(reduced_state::ReducedConstantTerm, basis::MPS, position::Int) + # Save the last `L` that is made to help with caching + # for DiskProjMPO + ll = reduced_state.lpos + if ll ≥ position + # Special case when nothing has to be done. + # Still need to change the position if lproj is + # being moved backward. + reduced_state.lpos = position + return nothing + end + # Make sure ll is at least 0 for the generic logic below + ll = max(ll, 0) + L = lproj(reduced_state) + while ll < position + L = L * basis[ll + 1] * dag(prime(reduced_state.state[ll + 1], "Link")) + reduced_state.environments[ll + 1] = L + ll += 1 + end + # Needed when moving lproj backward. + reduced_state.lpos = position + return reduced_state +end + +function ITensorMPS.makeR!(reduced_state::ReducedConstantTerm, basis::MPS, position::Int) + # Save the last `R` that is made to help with caching + # for DiskProjMPO + environment_position = reduced_state.rpos + if environment_position ≤ position + # Special case when nothing has to be done. + # Still need to change the position if rproj is + # being moved backward. + reduced_state.rpos = position + return nothing + end + N = length(reduced_state.state) + # Make sure environment_position is no bigger than `N + 1` for the generic logic below + environment_position = min(environment_position, N + 1) + right_environment = rproj(reduced_state) + while environment_position > position + right_environment = + right_environment * + basis[environment_position - 1] * + dag(prime(reduced_state.state[environment_position - 1], "Link")) + reduced_state.environments[environment_position - 1] = right_environment + environment_position -= 1 + end + reduced_state.rpos = position + return reduced_state +end + +function ITensors.contract(reduced_state::ReducedConstantTerm, v::ITensor) + reduced_state_tensors = Union{ITensor,OneITensor}[lproj(reduced_state)] + append!( + reduced_state_tensors, + [prime(t, "Link") for t in reduced_state.state[site_range(reduced_state)]], + ) + push!(reduced_state_tensors, rproj(reduced_state)) + + # 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(reduced_state_tensors)) == 1 + reverse!(reduced_state_tensors) + end + + # Apply the map + inner = v + for t in reduced_state_tensors + inner *= t + end + return inner +end + +# Contract the reduced constant term down to a since ITensor. +function ITensors.contract(reduced_state::ReducedConstantTerm) + reduced_state_tensors = Union{ITensor,OneITensor}[lproj(reduced_state)] + append!( + reduced_state_tensors, + [dag(prime(t, "Link")) for t in reduced_state.state[site_range(reduced_state)]], + ) + push!(reduced_state_tensors, rproj(reduced_state)) + + # 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(reduced_state_tensors)) == 1 + reverse!(reduced_state_tensors) + end + + # Apply the map + contracted_reduced_state = ITensor(true) + for t in reduced_state_tensors + contracted_reduced_state *= t + end + return contracted_reduced_state +end diff --git a/src/solvers/reducedcontractproblem.jl b/src/solvers/reducedcontractproblem.jl new file mode 100644 index 0000000..916d60c --- /dev/null +++ b/src/solvers/reducedcontractproblem.jl @@ -0,0 +1,122 @@ +using ITensors: ITensor + +""" +A ReducedContractProblem represents the application of an +MPO `operator` onto an MPS `input_state` but "projected" by +the basis of a different MPS `state` (which +could be an approximation to operator|state>). + +As an implementation of the AbstractProjMPO +type, it supports multiple `nsite` values for +one- and two-site algorithms. + +``` + *--*--*- -*--*--*--*--*--* +``` +""" +mutable struct ReducedContractProblem <: AbstractProjMPO + lpos::Int + rpos::Int + nsite::Int + input_state::MPS + operator::MPO + environments::Vector{ITensor} +end + +function ReducedContractProblem(input_state::MPS, operator::MPO) + lpos = 0 + rpos = length(operator) + 1 + nsite = 2 + environments = Vector{ITensor}(undef, length(operator)) + return ReducedContractProblem(lpos, rpos, nsite, input_state, operator, environments) +end + +function Base.getproperty(reduced_operator::ReducedContractProblem, sym::Symbol) + # This is for compatibility with `AbstractProjMPO`. + # TODO: Don't use `reduced_operator.H`, `reduced_operator.LR`, etc. + # in `AbstractProjMPO`. + if sym === :H + return getfield(reduced_operator, :operator) + elseif sym == :LR + return getfield(reduced_operator, :environments) + end + return getfield(reduced_operator, sym) +end + +function Base.copy(reduced_operator::ReducedContractProblem) + return ReducedContractProblem( + reduced_operator.lpos, + reduced_operator.rpos, + reduced_operator.nsite, + copy(reduced_operator.input_state), + copy(reduced_operator.operator), + copy(reduced_operator.environments), + ) +end + +Base.length(reduced_operator::ReducedContractProblem) = length(reduced_operator.operator) + +function ITensorMPS.set_nsite!(reduced_operator::ReducedContractProblem, nsite) + reduced_operator.nsite = nsite + return reduced_operator +end + +function ITensorMPS.makeL!(reduced_operator::ReducedContractProblem, state::MPS, k::Int) + # Save the last `L` that is made to help with caching + # for DiskProjMPO + ll = reduced_operator.lpos + if ll ≥ k + # Special case when nothing has to be done. + # Still need to change the position if lproj is + # being moved backward. + reduced_operator.lpos = k + return nothing + end + # Make sure ll is at least 0 for the generic logic below + ll = max(ll, 0) + L = lproj(reduced_operator) + while ll < k + L = + L * + reduced_operator.input_state[ll + 1] * + reduced_operator.operator[ll + 1] * + dag(state[ll + 1]) + reduced_operator.environments[ll + 1] = L + ll += 1 + end + # Needed when moving lproj backward. + reduced_operator.lpos = k + return reduced_operator +end + +function ITensorMPS.makeR!(reduced_operator::ReducedContractProblem, state::MPS, k::Int) + # Save the last `R` that is made to help with caching + # for DiskProjMPO + rl = reduced_operator.rpos + if rl ≤ k + # Special case when nothing has to be done. + # Still need to change the position if rproj is + # being moved backward. + reduced_operator.rpos = k + return nothing + end + N = length(reduced_operator.operator) + # Make sure rl is no bigger than `N + 1` for the generic logic below + rl = min(rl, N + 1) + R = rproj(reduced_operator) + while rl > k + R = + R * + reduced_operator.input_state[rl - 1] * + reduced_operator.operator[rl - 1] * + dag(state[rl - 1]) + reduced_operator.environments[rl - 1] = R + rl -= 1 + end + reduced_operator.rpos = k + return reduced_operator +end diff --git a/src/solvers/reducedlinearproblem.jl b/src/solvers/reducedlinearproblem.jl new file mode 100644 index 0000000..283fdfc --- /dev/null +++ b/src/solvers/reducedlinearproblem.jl @@ -0,0 +1,61 @@ +using ITensors: contract + +mutable struct ReducedLinearProblem <: AbstractProjMPO + reduced_operator::ProjMPO + reduced_constant_terms::Vector{ReducedConstantTerm} +end + +# Linear problem updater interface. +operator(reduced_problem::ReducedLinearProblem) = reduced_problem.reduced_operator +function constant_term(reduced_problem::ReducedLinearProblem) + constant_terms = map(reduced_problem.reduced_constant_terms) do reduced_constant_term + return contract(reduced_constant_term) + end + return dag(only(constant_terms)) +end + +function ReducedLinearProblem(operator::MPO, constant_term::MPS) + return ReducedLinearProblem(ProjMPO(operator), [ReducedConstantTerm(constant_term)]) +end + +function ReducedLinearProblem(operator::MPO, constant_terms::Vector{MPS}) + return ReducedLinearProblem(ProjMPO(operator), ReducedConstantTerm.(constant_terms)) +end + +function Base.copy(reduced_problem::ReducedLinearProblem) + return ReducedLinearProblem( + copy(reduced_problem.reduced_operator), copy(reduced_problem.reduced_constant_terms) + ) +end + +function ITensorMPS.nsite(reduced_problem::ReducedLinearProblem) + return nsite(reduced_problem.reduced_operator) +end + +function ITensorMPS.set_nsite!(reduced_problem::ReducedLinearProblem, nsite) + set_nsite!(reduced_problem.reduced_operator, nsite) + for m in reduced_problem.reduced_constant_terms + set_nsite!(m, nsite) + end + return reduced_problem +end + +function ITensorMPS.makeL!(reduced_problem::ReducedLinearProblem, state::MPS, position::Int) + makeL!(reduced_problem.reduced_operator, state, position) + for reduced_constant_term in reduced_problem.reduced_constant_terms + makeL!(reduced_constant_term, state, position) + end + return reduced_problem +end + +function ITensorMPS.makeR!(reduced_problem::ReducedLinearProblem, state::MPS, position::Int) + makeR!(reduced_problem.reduced_operator, state, position) + for reduced_constant_term in reduced_problem.reduced_constant_terms + makeR!(reduced_constant_term, state, position) + end + return reduced_problem +end + +function ITensors.contract(reduced_problem::ReducedLinearProblem, v::ITensor) + return contract(reduced_problem.reduced_operator, v) +end diff --git a/src/solvers/sweep_update.jl b/src/solvers/sweep_update.jl new file mode 100644 index 0000000..46ffa3e --- /dev/null +++ b/src/solvers/sweep_update.jl @@ -0,0 +1,476 @@ +using ITensors: ITensors, uniqueinds +using LinearAlgebra: norm, normalize!, svd +using Printf: @printf + +function sweep_update( + order::TDVPOrder, + reduced_operator, + state::MPS; + current_time=nothing, + time_step=nothing, + kwargs..., +) + order_orderings = orderings(order) + order_sub_time_steps = sub_time_steps(order) + if !isnothing(time_step) + order_sub_time_steps = eltype(time_step).(order_sub_time_steps) + order_sub_time_steps *= time_step + end + info = nothing + sub_time_step = nothing + for substep in 1:length(order_sub_time_steps) + if !isnothing(time_step) + sub_time_step = order_sub_time_steps[substep] + end + state, reduced_operator, info = sub_sweep_update( + order_orderings[substep], + reduced_operator, + state; + current_time, + time_step=sub_time_step, + kwargs..., + ) + if !isnothing(time_step) + current_time += sub_time_step + end + end + return state, reduced_operator, info +end + +isforward(direction::Base.ForwardOrdering) = true +isforward(direction::Base.ReverseOrdering) = false +isreverse(direction) = !isforward(direction) + +function sweep_bonds(direction::Base.ForwardOrdering, n::Int; ncenter::Int) + return 1:(n - ncenter + 1) +end + +function sweep_bonds(direction::Base.ReverseOrdering, n::Int; ncenter::Int) + return reverse(sweep_bonds(Base.Forward, n; ncenter)) +end + +is_forward_done(direction::Base.ForwardOrdering, b, n; ncenter) = (b + ncenter - 1 == n) +is_forward_done(direction::Base.ReverseOrdering, b, n; ncenter) = false +is_reverse_done(direction::Base.ForwardOrdering, b, n; ncenter) = false +is_reverse_done(direction::Base.ReverseOrdering, b, n; ncenter) = (b == 1) +function is_half_sweep_done(direction, b, n; ncenter) + return is_forward_done(direction, b, n; ncenter) || + is_reverse_done(direction, b, n; ncenter) +end + +function sub_sweep_update( + direction::Base.Ordering, + reduced_operator, + state::MPS; + updater, + updater_kwargs, + which_decomp=nothing, + svd_alg=nothing, + sweep=default_sweep(), + current_time=nothing, + time_step=nothing, + nsite=default_nsite(), + reverse_step=default_reverse_step(), + normalize=default_normalize(), + (observer!)=default_observer(), + outputlevel=default_outputlevel(), + maxdim=default_maxdim(), + mindim=default_mindim(), + cutoff=default_cutoff(ITensors.scalartype(state)), + noise=default_noise(), +) + reduced_operator = copy(reduced_operator) + state = copy(state) + if length(state) == 1 + error( + "`tdvp`, `dmrg`, `linsolve`, etc. currently does not support system sizes of 1. You can diagonalize the MPO tensor directly with tools like `LinearAlgebra.eigen`, `KrylovKit.exponentiate`, etc.", + ) + end + N = length(state) + set_nsite!(reduced_operator, nsite) + if isforward(direction) + if !isortho(state) || orthocenter(state) != 1 + orthogonalize!(state, 1) + end + @assert isortho(state) && orthocenter(state) == 1 + position!(reduced_operator, state, 1) + elseif isreverse(direction) + if !isortho(state) || orthocenter(state) != N - nsite + 1 + orthogonalize!(state, N - nsite + 1) + end + @assert(isortho(state) && (orthocenter(state) == N - nsite + 1)) + position!(reduced_operator, state, N - nsite + 1) + end + maxtruncerr = 0.0 + info = nothing + for b in sweep_bonds(direction, N; ncenter=nsite) + current_time, maxtruncerr, spec, info = region_update!( + reduced_operator, + state, + b; + updater, + updater_kwargs, + nsite, + reverse_step, + current_time, + outputlevel, + time_step, + normalize, + direction, + noise, + which_decomp, + svd_alg, + cutoff, + maxdim, + mindim, + maxtruncerr, + ) + if outputlevel >= 2 + if nsite == 1 + @printf("Sweep %d, direction %s, bond (%d,) \n", sweep, direction, b) + elseif nsite == 2 + @printf("Sweep %d, direction %s, bond (%d,%d) \n", sweep, direction, b, b + 1) + end + print(" Truncated using") + @printf(" cutoff=%.1E", cutoff) + @printf(" maxdim=%.1E", maxdim) + print(" mindim=", mindim) + print(" current_time=", round(current_time; digits=3)) + println() + if spec != nothing + @printf( + " Trunc. err=%.2E, bond dimension %d\n", spec.truncerr, dim(linkind(state, b)) + ) + end + flush(stdout) + end + update_observer!( + observer!; + state, + reduced_operator, + bond=b, + sweep, + half_sweep=isforward(direction) ? 1 : 2, + spec, + outputlevel, + half_sweep_is_done=is_half_sweep_done(direction, b, N; ncenter=nsite), + current_time, + info, + ) + end + # Just to be sure: + normalize && normalize!(state) + return state, reduced_operator, (; maxtruncerr) +end + +function region_update!( + reduced_operator, + state, + b; + updater, + updater_kwargs, + nsite, + reverse_step, + current_time, + outputlevel, + time_step, + normalize, + direction, + noise, + which_decomp, + svd_alg, + cutoff, + maxdim, + mindim, + maxtruncerr, +) + return region_update!( + Val(nsite), + Val(reverse_step), + reduced_operator, + state, + b; + updater, + updater_kwargs, + current_time, + outputlevel, + time_step, + normalize, + direction, + noise, + which_decomp, + svd_alg, + cutoff, + maxdim, + mindim, + maxtruncerr, + ) +end + +function region_update!( + nsite_val::Val{1}, + reverse_step_val::Val{false}, + reduced_operator, + state, + b; + updater, + updater_kwargs, + current_time, + outputlevel, + time_step, + normalize, + direction, + noise, + which_decomp, + svd_alg, + cutoff, + maxdim, + mindim, + maxtruncerr, +) + N = length(state) + nsite = 1 + # Do 'forwards' evolution step + set_nsite!(reduced_operator, nsite) + position!(reduced_operator, state, b) + reduced_state = state[b] + internal_kwargs = (; current_time, time_step, outputlevel) + reduced_state, info = updater( + reduced_operator, reduced_state; internal_kwargs, updater_kwargs... + ) + if !isnothing(time_step) + current_time += time_step + end + normalize && (reduced_state /= norm(reduced_state)) + spec = nothing + state[b] = reduced_state + if !is_half_sweep_done(direction, b, N; ncenter=nsite) + # Move ortho center + Δ = (isforward(direction) ? +1 : -1) + orthogonalize!(state, b + Δ) + end + return current_time, maxtruncerr, spec, info +end + +function region_update!( + nsite_val::Val{1}, + reverse_step_val::Val{true}, + reduced_operator, + state, + b; + updater, + updater_kwargs, + current_time, + outputlevel, + time_step, + normalize, + direction, + noise, + which_decomp, + svd_alg, + cutoff, + maxdim, + mindim, + maxtruncerr, +) + N = length(state) + nsite = 1 + # Do 'forwards' evolution step + set_nsite!(reduced_operator, nsite) + position!(reduced_operator, state, b) + reduced_state = state[b] + internal_kwargs = (; current_time, time_step, outputlevel) + reduced_state, info = updater( + reduced_operator, reduced_state; internal_kwargs, updater_kwargs... + ) + current_time += time_step + normalize && (reduced_state /= norm(reduced_state)) + spec = nothing + state[b] = reduced_state + if !is_half_sweep_done(direction, b, N; ncenter=nsite) + # Do backwards evolution step + b1 = (isforward(direction) ? b + 1 : b) + Δ = (isforward(direction) ? +1 : -1) + uinds = uniqueinds(reduced_state, state[b + Δ]) + U, S, V = svd(reduced_state, uinds) + state[b] = U + bond_reduced_state = S * V + if isforward(direction) + ITensorMPS.setleftlim!(state, b) + elseif isreverse(direction) + ITensorMPS.setrightlim!(state, b) + end + set_nsite!(reduced_operator, nsite - 1) + position!(reduced_operator, state, b1) + internal_kwargs = (; current_time, time_step=-time_step, outputlevel) + bond_reduced_state, info = updater( + reduced_operator, bond_reduced_state; internal_kwargs, updater_kwargs... + ) + current_time -= time_step + normalize && (bond_reduced_state ./= norm(bond_reduced_state)) + state[b + Δ] = bond_reduced_state * state[b + Δ] + if isforward(direction) + ITensorMPS.setrightlim!(state, b + Δ + 1) + elseif isreverse(direction) + ITensorMPS.setleftlim!(state, b + Δ - 1) + end + set_nsite!(reduced_operator, nsite) + end + return current_time, maxtruncerr, spec, info +end + +function region_update!( + nsite_val::Val{2}, + reverse_step_val::Val{false}, + reduced_operator, + state, + b; + updater, + updater_kwargs, + current_time, + time_step, + outputlevel, + normalize, + direction, + noise, + which_decomp, + svd_alg, + cutoff, + maxdim, + mindim, + maxtruncerr, +) + N = length(state) + nsite = 2 + # Do 'forwards' evolution step + set_nsite!(reduced_operator, nsite) + position!(reduced_operator, state, b) + reduced_state = state[b] * state[b + 1] + internal_kwargs = (; current_time, time_step, outputlevel) + reduced_state, info = updater( + reduced_operator, reduced_state; internal_kwargs, updater_kwargs... + ) + if !isnothing(time_step) + current_time += time_step + end + normalize && (reduced_state /= norm(reduced_state)) + spec = nothing + ortho = isforward(direction) ? "left" : "right" + drho = nothing + if noise > 0.0 && isforward(direction) + drho = noise * noiseterm(reduced_operator, reduced_state, ortho) + end + spec = replacebond!( + state, + b, + reduced_state; + maxdim, + mindim, + cutoff, + eigen_perturbation=drho, + ortho=ortho, + normalize, + which_decomp, + svd_alg, + ) + maxtruncerr = max(maxtruncerr, spec.truncerr) + return current_time, maxtruncerr, spec, info +end + +function region_update!( + nsite_val::Val{2}, + reverse_step_val::Val{true}, + reduced_operator, + state, + b; + updater, + updater_kwargs, + current_time, + time_step, + outputlevel, + normalize, + direction, + noise, + which_decomp, + svd_alg, + cutoff, + maxdim, + mindim, + maxtruncerr, +) + N = length(state) + nsite = 2 + # Do 'forwards' evolution step + set_nsite!(reduced_operator, nsite) + position!(reduced_operator, state, b) + reduced_state = state[b] * state[b + 1] + internal_kwargs = (; current_time, time_step, outputlevel) + reduced_state, info = updater( + reduced_operator, reduced_state; internal_kwargs, updater_kwargs... + ) + current_time += time_step + normalize && (reduced_state /= norm(reduced_state)) + spec = nothing + ortho = isforward(direction) ? "left" : "right" + drho = nothing + if noise > 0.0 && isforward(direction) + drho = noise * noiseterm(reduced_operator, phi, ortho) + end + spec = replacebond!( + state, + b, + reduced_state; + maxdim, + mindim, + cutoff, + eigen_perturbation=drho, + ortho=ortho, + normalize, + which_decomp, + svd_alg, + ) + maxtruncerr = max(maxtruncerr, spec.truncerr) + if !is_half_sweep_done(direction, b, N; ncenter=nsite) + # Do backwards evolution step + b1 = (isforward(direction) ? b + 1 : b) + Δ = (isforward(direction) ? +1 : -1) + bond_reduced_state = state[b1] + set_nsite!(reduced_operator, nsite - 1) + position!(reduced_operator, state, b1) + internal_kwargs = (; current_time, time_step=-time_step, outputlevel) + bond_reduced_state, info = updater( + reduced_operator, bond_reduced_state; internal_kwargs, updater_kwargs... + ) + current_time -= time_step + normalize && (bond_reduced_state /= norm(bond_reduced_state)) + state[b1] = bond_reduced_state + set_nsite!(reduced_operator, nsite) + end + return current_time, maxtruncerr, spec, info +end + +function region_!( + ::Val{nsite}, + ::Val{reverse_step}, + reduced_operator, + state, + b; + updater, + updater_kwargs, + current_time, + outputlevel, + time_step, + normalize, + direction, + noise, + which_decomp, + svd_alg, + cutoff, + maxdim, + mindim, + maxtruncerr, +) where {nsite,reverse_step} + return error( + "`tdvp`, `dmrg`, `linsolve`, etc. with `nsite=$nsite` and `reverse_step=$reverse_step` not implemented.", + ) +end diff --git a/src/solvers/tdvp.jl b/src/solvers/tdvp.jl new file mode 100644 index 0000000..651bd34 --- /dev/null +++ b/src/solvers/tdvp.jl @@ -0,0 +1,92 @@ +using ITensors: Algorithm, @Algorithm_str +using KrylovKit: exponentiate + +function exponentiate_updater(operator, init; internal_kwargs, kwargs...) + state, info = exponentiate(operator, internal_kwargs.time_step, init; kwargs...) + return state, (; info) +end + +function applyexp_updater(operator, init; internal_kwargs, kwargs...) + state, info = applyexp(operator, internal_kwargs.time_step, init; kwargs...) + return state, (; info) +end + +tdvp_updater(updater_backend::String) = tdvp_updater(Algorithm(updater_backend)) +tdvp_updater(::Algorithm"exponentiate") = exponentiate_updater +tdvp_updater(::Algorithm"applyexp") = applyexp_updater +function tdvp_updater(updater_backend::Algorithm) + return error("`updater_backend=$(String(updater_backend))` not recognized.") +end + +function time_step_and_nsteps(t, time_step::Nothing, nsteps::Nothing) + # Default to 1 step. + nsteps = 1 + return time_step_and_nsteps(t, time_step, nsteps) +end + +function time_step_and_nsteps(t, time_step::Nothing, nsteps) + return t / nsteps, nsteps +end + +function time_step_and_nsteps(t, time_step, nsteps::Nothing) + nsteps_float = t / time_step + nsteps_rounded = round(nsteps_float) + if abs(nsteps_float - nsteps_rounded) ≉ 0 + return error("`t / time_step = $t / $time_step = $(t / time_step)` must be an integer.") + end + return time_step, Int(nsteps_rounded) +end + +function time_step_and_nsteps(t, time_step, nsteps) + if time_step * nsteps ≠ t + return error( + "Calling `tdvp(operator, t, state; time_step, nsteps, kwargs...)` with `t = $t`, `time_step = $time_step`, and `nsteps = $nsteps` must satisfy `time_step * nsteps == t`, while `time_step * nsteps = $time_step * $nsteps = $(time_step * nsteps)`.", + ) + end + return time_step, nsteps +end + +""" + tdvp(operator, t::Number, init::MPS; time_step, nsteps, kwargs...) + +Use the time dependent variational principle (TDVP) algorithm +to compute `exp(t * operator) * init` using an efficient algorithm based +on alternating optimization of the MPS tensors and local Krylov +exponentiation of `operator`. + +Specify one of `time_step` or `nsteps`. If they are both specified, they +must satisfy `time_step * nsteps == t`. If neither are specified, the +default is `nsteps=1`, which means that `time_step == t`. + +Returns: +* `state::MPS` - time-evolved MPS + +""" +function tdvp( + operator, + t::Number, + init::MPS; + updater_backend="exponentiate", + updater=tdvp_updater(updater_backend), + reverse_step=true, + time_step=nothing, + time_start=zero(t), + nsweeps=nothing, + nsteps=nsweeps, + (step_observer!)=default_sweep_observer(), + (sweep_observer!)=step_observer!, + kwargs..., +) + time_step, nsteps = time_step_and_nsteps(t, time_step, nsteps) + return alternating_update( + operator, + init; + updater, + reverse_step, + nsweeps=nsteps, + time_start, + time_step, + sweep_observer!, + kwargs..., + ) +end diff --git a/src/solvers/tdvporder.jl b/src/solvers/tdvporder.jl new file mode 100644 index 0000000..e37e00e --- /dev/null +++ b/src/solvers/tdvporder.jl @@ -0,0 +1,24 @@ +struct TDVPOrder{order,direction} end + +TDVPOrder(order::Int, direction::Base.Ordering) = TDVPOrder{order,direction}() + +orderings(::TDVPOrder) = error("Not implemented") +sub_time_steps(::TDVPOrder) = error("Not implemented") + +function orderings(::TDVPOrder{1,direction}) where {direction} + return [direction, Base.ReverseOrdering(direction)] +end +sub_time_steps(::TDVPOrder{1}) = [1, 0] + +function orderings(::TDVPOrder{2,direction}) where {direction} + return [direction, Base.ReverseOrdering(direction)] +end +sub_time_steps(::TDVPOrder{2}) = [1 / 2, 1 / 2] + +function orderings(::TDVPOrder{4,direction}) where {direction} + return [direction, Base.ReverseOrdering(direction)] +end +function sub_time_steps(::TDVPOrder{4}) + s = 1 / (2 - 2^(1 / 3)) + return [s / 2, s / 2, (1 - 2 * s) / 2, (1 - 2 * s) / 2, s / 2, s / 2] +end diff --git a/src/solvers/timedependentsum.jl b/src/solvers/timedependentsum.jl new file mode 100644 index 0000000..f2998de --- /dev/null +++ b/src/solvers/timedependentsum.jl @@ -0,0 +1,68 @@ +using ITensors: ITensor, inds, permute + +# Represents a time-dependent sum of terms: +# +# expr(t) = coefficients(expr)[1](t) * terms(expr)[1] + coefficients(expr)[2](t) * terms(expr)[2] + … +# +struct TimeDependentSum{Coefficients,Terms} + coefficients::Coefficients + terms::Terms +end + +coefficients(expr::TimeDependentSum) = expr.coefficients +terms(expr::TimeDependentSum) = expr.terms +function Base.copy(expr::TimeDependentSum) + return TimeDependentSum(coefficients(expr), copy.(terms(expr))) +end + +function Base.:*(c::Number, expr::TimeDependentSum) + scaled_coefficients = map(coefficient -> (t -> c * coefficient(t)), coefficients(expr)) + return TimeDependentSum(scaled_coefficients, terms(expr)) +end +Base.:*(expr::TimeDependentSum, c::Number) = c * expr + +# Evaluating a `TimeDependentSum` at a certain time +# returns a `ScaledSum` at that time. +function (expr::TimeDependentSum)(t::Number) + coefficients_t = map(coefficient -> coefficient(t), coefficients(expr)) + return ScaledSum(coefficients_t, terms(expr)) +end + +# alternating_update inteface +function reduced_operator(operator::TimeDependentSum) + return TimeDependentSum(coefficients(operator), reduced_operator.(terms(operator))) +end +function ITensorMPS.set_nsite!(operator::TimeDependentSum, nsite) + foreach(t -> set_nsite!(t, nsite), terms(operator)) + return operator +end +function ITensorMPS.position!(operator::TimeDependentSum, state, position) + foreach(t -> position!(t, state, position), terms(operator)) + return operator +end + +# Represents the sum of scaled terms: +# +# H = coefficient[1] * H[1] + coefficient * H[2] + … +# +struct ScaledSum{Coefficients,Terms} + coefficients::Coefficients + terms::Terms +end + +coefficients(expr::ScaledSum) = expr.coefficients +terms(expr::ScaledSum) = expr.terms + +# Apply the scaled sum of terms: +# +# expr(x) = coefficients(expr)[1] * terms(expr)[1](x) + coefficients(expr)[2] * terms(expr)[2](x) + … +# +# onto x. +function scaledsum_apply(expr, x) + return mapreduce(+, zip(coefficients(expr), terms(expr))) do coefficient_and_term + coefficient, term = coefficient_and_term + return coefficient * term(x) + end +end +(expr::ScaledSum)(x) = scaledsum_apply(expr, x) +(expr::ScaledSum)(x::ITensor) = permute(scaledsum_apply(expr, x), inds(x)) diff --git a/src/solvers/update_observer.jl b/src/solvers/update_observer.jl new file mode 100644 index 0000000..3d388ab --- /dev/null +++ b/src/solvers/update_observer.jl @@ -0,0 +1,24 @@ +struct EmptyObserver end +update_observer!(observer::EmptyObserver; kwargs...) = observer + +struct ValuesObserver{Values<:NamedTuple} + values::Values +end +function update_observer!(observer::ValuesObserver; kwargs...) + for key in keys(observer.values) + observer.values[key][] = kwargs[key] + end + return observer +end +values_observer(; kwargs...) = ValuesObserver(NamedTuple(kwargs)) + +struct ComposedObservers{Observers<:Tuple} + observers::Observers +end +compose_observers(observers...) = ComposedObservers(observers) +function update_observer!(observer::ComposedObservers; kwargs...) + for observerᵢ in observer.observers + update_observer!(observerᵢ; kwargs...) + end + return observer +end diff --git a/test/Ops/Project.toml b/test/Ops/Project.toml new file mode 100644 index 0000000..a355267 --- /dev/null +++ b/test/Ops/Project.toml @@ -0,0 +1,3 @@ +[deps] +ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" +ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" diff --git a/test/Ops/test_ops_mpo.jl b/test/Ops/test_ops_mpo.jl index 1348aae..c04faea 100644 --- a/test/Ops/test_ops_mpo.jl +++ b/test/Ops/test_ops_mpo.jl @@ -1,9 +1,9 @@ @eval module $(gensym()) -using Test -using ITensors -using ITensors.Ops -using ITensors.ITensorMPS: ITensorMPS -using LinearAlgebra +using ITensors: ITensor, contract, op, replaceprime +using ITensors.Ops: Op, OpSum, Prod, Scaled, Sum, expand +using ITensorMPS: ITensorMPS, MPO, siteinds +using LinearAlgebra: I, norm +using Test: @test, @testset @testset "Ops to MPO" begin ∑H = Sum{Op}() diff --git a/test/Ops/test_trotter.jl b/test/Ops/test_trotter.jl index 54554df..f03ea7d 100644 --- a/test/Ops/test_trotter.jl +++ b/test/Ops/test_trotter.jl @@ -1,7 +1,8 @@ @eval module $(gensym()) -using Test -using ITensors -using ITensors.Ops +using Test: @test, @testset +using ITensorMPS: MPO, MPS, siteinds +using ITensors: ITensor, apply, contract, replaceprime +using ITensors.Ops: Op, Prod, Sum, Trotter function heisenberg(N) os = Sum{Op}() diff --git a/test/backup/test_exports.jl b/test/base/test_exports.jl similarity index 100% rename from test/backup/test_exports.jl rename to test/base/test_exports.jl diff --git a/test/backup/Project.toml b/test/base/test_solvers/Project.toml similarity index 80% rename from test/backup/Project.toml rename to test/base/test_solvers/Project.toml index 824fb9a..dcd1d96 100644 --- a/test/backup/Project.toml +++ b/test/base/test_solvers/Project.toml @@ -1,12 +1,16 @@ [deps] Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" -ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensorTDVP = "25707e16-a4db-4a07-99d9-4d67b7af0342" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0" OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[compat] +ITensors = "0.6.7" diff --git a/test/backup/runtests.jl b/test/base/test_solvers/runtests.jl similarity index 72% rename from test/backup/runtests.jl rename to test/base/test_solvers/runtests.jl index 44e961e..5691f1f 100644 --- a/test/backup/runtests.jl +++ b/test/base/test_solvers/runtests.jl @@ -1,11 +1,11 @@ @eval module $(gensym()) using Test: @testset -using ITensorMPS: ITensorMPS -test_path = joinpath(pkgdir(ITensorMPS), "test") +using ITensorTDVP: ITensorTDVP +test_path = joinpath(pkgdir(ITensorTDVP), "test") test_files = filter( file -> startswith(file, "test_") && endswith(file, ".jl"), readdir(test_path) ) -@testset "ITensorMPS.jl" begin +@testset "ITensorTDVP.jl" begin @testset "$filename" for filename in test_files println("Running $filename") @time include(joinpath(test_path, filename)) diff --git a/test/base/test_solvers/test_contract.jl b/test/base/test_solvers/test_contract.jl new file mode 100644 index 0000000..e0dc288 --- /dev/null +++ b/test/base/test_solvers/test_contract.jl @@ -0,0 +1,66 @@ +@eval module $(gensym()) +using ITensors: ITensors, dag, delta, denseblocks +using ITensors: MPO, OpSum, apply, contract, inner, random_mps, siteinds, truncate! +using ITensorTDVP: ITensorTDVP +using StableRNGs: StableRNG +using Test: @test, @test_throws, @testset +@testset "Contract MPO (eltype=$elt, conserve_qns=$conserve_qns)" for elt in ( + Float32, Float64, Complex{Float32}, Complex{Float64} + ), + conserve_qns in [false, true] + + N = 20 + s = siteinds("S=1/2", N; conserve_qns) + rng = StableRNG(1234) + psi = random_mps(rng, elt, s, j -> isodd(j) ? "↑" : "↓"; linkdims=8) + os = OpSum() + for j in 1:(N - 1) + os += 0.5, "S+", j, "S-", j + 1 + os += 0.5, "S-", j, "S+", j + 1 + os += "Sz", j, "Sz", j + 1 + end + for j in 1:(N - 2) + os += 0.5, "S+", j, "S-", j + 2 + os += 0.5, "S-", j, "S+", j + 2 + os += "Sz", j, "Sz", j + 2 + end + H = MPO(elt, os, s) + @testset "apply (standard indices, nsite=2)" begin + Hpsi = apply(H, psi; alg="fit", nsweeps=2) + @test_throws ErrorException apply(H, psi; alg="fit") + @test ITensors.scalartype(Hpsi) == elt + @test inner(psi, Hpsi) ≈ inner(psi', H, psi) rtol = 10 * √eps(real(elt)) + end + @testset "contract (non-standard indices)" begin + # Change "top" indices of MPO to be a different set + t = siteinds("S=1/2", N; conserve_qns) + Ht = deepcopy(H) + psit = deepcopy(psi) + for j in 1:N + Ht[j] *= delta(elt, dag(s[j])', t[j]) + psit[j] *= delta(elt, dag(s[j]), t[j]) + end + + # Test with nsweeps=2 + Hpsit = contract(Ht, psi; alg="fit", nsweeps=2) + @test ITensors.scalartype(Hpsit) == elt + @test inner(psit, Hpsit) ≈ inner(psit, Ht, psi) rtol = 10 * √eps(real(elt)) + + # Test with less good initial guess MPS not equal to psi + psit_guess = copy(psit) + truncate!(psit_guess; maxdim=2) + Hpsit = contract(Ht, psi; alg="fit", nsweeps=4, init=psit_guess) + @test ITensors.scalartype(Hpsit) == elt + @test inner(psit, Hpsit) ≈ inner(psit, Ht, psi) rtol = 20 * √eps(real(elt)) + end + @testset "apply (standard indices, nsite=1)" begin + # Test with nsite=1 + Hpsi_guess = apply(H, psi; alg="naive", cutoff=1e-4) + Hpsi = apply(H, psi; alg="fit", init=Hpsi_guess, nsite=1, nsweeps=2) + @test ITensors.scalartype(Hpsi) == elt + scale(::Type{Float32}) = 10^2 + scale(::Type{Float64}) = 10^6 + @test inner(psi, Hpsi) ≈ inner(psi', H, psi) rtol = √eps(real(elt)) * scale(real(elt)) + end +end +end diff --git a/test/base/test_solvers/test_dmrg.jl b/test/base/test_solvers/test_dmrg.jl new file mode 100644 index 0000000..a7237c8 --- /dev/null +++ b/test/base/test_solvers/test_dmrg.jl @@ -0,0 +1,37 @@ +@eval module $(gensym()) +using ITensors: ITensors, MPO, OpSum, inner, random_mps, siteinds +using ITensorTDVP: ITensorTDVP +using StableRNGs: StableRNG +using Test: @test, @test_throws, @testset +@testset "DMRG (eltype=$elt, nsite=$nsite, conserve_qns=$conserve_qns)" for elt in ( + Float32, Float64, Complex{Float32}, Complex{Float64} + ), + nsite in [1, 2], + conserve_qns in [false, true] + + N = 10 + cutoff = eps(real(elt)) * 10^4 + s = siteinds("S=1/2", N; conserve_qns) + os = OpSum() + for j in 1:(N - 1) + os += 0.5, "S+", j, "S-", j + 1 + os += 0.5, "S-", j, "S+", j + 1 + os += "Sz", j, "Sz", j + 1 + end + H = MPO(elt, os, s) + rng = StableRNG(1234) + psi = random_mps(rng, elt, s, j -> isodd(j) ? "↑" : "↓"; linkdims=20) + nsweeps = 10 + maxdim = [10, 20, 40, 100] + @test_throws ErrorException ITensorTDVP.dmrg(H, psi; maxdim, cutoff, nsite) + e, psi = ITensorTDVP.dmrg( + H, psi; nsweeps, maxdim, cutoff, nsite, updater_kwargs=(; krylovdim=3, maxiter=1) + ) + @test inner(psi', H, psi) ≈ e + e2, psi2 = ITensors.dmrg(H, psi; nsweeps, maxdim, cutoff, outputlevel=0) + @test ITensors.scalartype(psi2) == elt + @test e2 isa real(elt) + @test e ≈ e2 rtol = √(eps(real(elt))) * 10 + @test inner(psi', H, psi) ≈ inner(psi2', H, psi2) rtol = √(eps(real(elt))) * 10 +end +end diff --git a/test/base/test_solvers/test_dmrg_x.jl b/test/base/test_solvers/test_dmrg_x.jl new file mode 100644 index 0000000..17c8cf3 --- /dev/null +++ b/test/base/test_solvers/test_dmrg_x.jl @@ -0,0 +1,55 @@ +@eval module $(gensym()) +using ITensors: ITensors, MPO, MPS, OpSum, ProjMPO, inner, siteinds +using ITensorTDVP: dmrg_x +using Random: Random +using StableRNGs: StableRNG +using Test: @test, @test_throws, @testset +@testset "DMRG-X (eltype=$elt, conserve_qns=$conserve_qns)" for elt in ( + Float32, Float64, Complex{Float32}, Complex{Float64} + ), + conserve_qns in [false, true] + + function heisenberg(n; h=zeros(n)) + os = OpSum() + for j in 1:(n - 1) + os += 0.5, "S+", j, "S-", j + 1 + os += 0.5, "S-", j, "S+", j + 1 + os += "Sz", j, "Sz", j + 1 + end + for j in 1:n + if h[j] ≠ 0 + os -= h[j], "Sz", j + end + end + return os + end + n = 10 + s = siteinds("S=1/2", n; conserve_qns) + Random.seed!(12) + W = 12 + # Random fields h ∈ [-W, W] + rng = StableRNG(1234) + h = W * (2 * rand(rng, real(elt), n) .- 1) + H = MPO(elt, heisenberg(n; h), s) + initstate = rand(rng, ["↑", "↓"], n) + ψ = MPS(elt, s, initstate) + @test_throws ErrorException dmrg_x(H, ψ; nsite=2, maxdim=20, cutoff=1e-10) + dmrg_x_kwargs = (; nsweeps=20, normalize=true, maxdim=20, cutoff=1e-10, outputlevel=0) + e, ϕ = dmrg_x(H, ψ; nsite=2, dmrg_x_kwargs...) + @test ITensors.scalartype(ϕ) == elt + @test inner(ϕ', H, ϕ) / inner(ϕ, ϕ) ≈ e + @test inner(H, ψ, H, ψ) ≉ inner(ψ', H, ψ)^2 rtol = √eps(real(elt)) + @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ', H, ϕ) / inner(ϕ, ϕ) rtol = 1e-1 + @test inner(H, ϕ, H, ϕ) ≈ inner(ϕ', H, ϕ)^2 rtol = √eps(real(elt)) + ẽ, ϕ̃ = dmrg_x(ProjMPO(H), ϕ; nsite=1, dmrg_x_kwargs...) + @test ITensors.scalartype(ϕ̃) == elt + @test inner(ϕ̃', H, ϕ̃) / inner(ϕ̃, ϕ̃) ≈ ẽ + @test inner(ψ', H, ψ) / inner(ψ, ψ) ≈ inner(ϕ̃', H, ϕ̃) / inner(ϕ̃, ϕ̃) rtol = 1e-1 + scale(::Type{Float32}) = 10^2 + scale(::Type{Float64}) = 10^5 + @test inner(H, ϕ̃, H, ϕ̃) ≈ inner(ϕ̃', H, ϕ̃)^2 rtol = + √(eps(real(elt))) * scale(real(elt)) + # Sometimes broken, sometimes not + # @test abs(loginner(ϕ̃, ϕ) / n) ≈ 0.0 atol = 1e-6 +end +end diff --git a/test/backup/test_examples.jl b/test/base/test_solvers/test_examples.jl similarity index 80% rename from test/backup/test_examples.jl rename to test/base/test_solvers/test_examples.jl index 3e0f481..88f0ac7 100644 --- a/test/backup/test_examples.jl +++ b/test/base/test_solvers/test_examples.jl @@ -1,12 +1,12 @@ @eval module $(gensym()) using Suppressor: @suppress -using ITensorMPS: ITensorMPS +using ITensorTDVP: ITensorTDVP using Test: @testset @testset "Run examples" begin examples_files = [ "01_tdvp.jl", "02_dmrg-x.jl", "03_tdvp_time_dependent.jl", "04_tdvp_observers.jl" ] - examples_path = joinpath(pkgdir(ITensorMPS), "examples") + examples_path = joinpath(pkgdir(ITensorTDVP), "examples") @testset "Running example file $f" for f in examples_files println("Running example file $f") @suppress include(joinpath(examples_path, f)) diff --git a/test/base/test_solvers/test_expand.jl b/test/base/test_solvers/test_expand.jl new file mode 100644 index 0000000..81b9d9a --- /dev/null +++ b/test/base/test_solvers/test_expand.jl @@ -0,0 +1,95 @@ +@eval module $(gensym()) +using ITensors: scalartype +using ITensors.ITensorMPS: + OpSum, MPO, MPS, inner, linkdims, maxlinkdim, random_mps, siteinds +using ITensorTDVP: dmrg, expand, tdvp +using LinearAlgebra: normalize +using StableRNGs: StableRNG +using Test: @test, @testset +const elts = (Float32, Float64, Complex{Float32}, Complex{Float64}) +@testset "expand (eltype=$elt)" for elt in elts + @testset "expand (alg=\"orthogonalize\", conserve_qns=$conserve_qns, eltype=$elt)" for conserve_qns in + ( + false, true + ) + n = 6 + s = siteinds("S=1/2", n; conserve_qns) + rng = StableRNG(1234) + state = random_mps(rng, elt, s, j -> isodd(j) ? "↑" : "↓"; linkdims=4) + reference = random_mps(rng, elt, s, j -> isodd(j) ? "↑" : "↓"; linkdims=2) + state_expanded = expand(state, [reference]; alg="orthogonalize") + @test scalartype(state_expanded) === elt + @test inner(state_expanded, state) ≈ inner(state, state) + @test inner(state_expanded, reference) ≈ inner(state, reference) + end + @testset "expand (alg=\"global_krylov\", conserve_qns=$conserve_qns, eltype=$elt)" for conserve_qns in + ( + false, true + ) + n = 10 + s = siteinds("S=1/2", n; conserve_qns) + opsum = OpSum() + for j in 1:(n - 1) + opsum += 0.5, "S+", j, "S-", j + 1 + opsum += 0.5, "S-", j, "S+", j + 1 + opsum += "Sz", j, "Sz", j + 1 + end + operator = MPO(elt, opsum, s) + state = MPS(elt, s, j -> isodd(j) ? "↑" : "↓") + state_expanded = expand(state, operator; alg="global_krylov") + @test scalartype(state_expanded) === elt + @test maxlinkdim(state_expanded) > 1 + @test inner(state_expanded, state) ≈ inner(state, state) + end + @testset "Decoupled ladder (alg=\"global_krylov\", eltype=$elt)" begin + nx = 10 + ny = 2 + n = nx * ny + s = siteinds("S=1/2", n) + opsum = OpSum() + for j in 1:2:(n - 2) + opsum += 1 / 2, "S+", j, "S-", j + 2 + opsum += 1 / 2, "S-", j, "S+", j + 2 + opsum += "Sz", j, "Sz", j + 2 + end + for j in 2:2:(n - 2) + opsum += 1 / 2, "S+", j, "S-", j + 2 + opsum += 1 / 2, "S-", j, "S+", j + 2 + opsum += "Sz", j, "Sz", j + 2 + end + operator = MPO(elt, opsum, s) + rng = StableRNG(1234) + init = random_mps(rng, elt, s; linkdims=30) + reference_energy, reference_state = dmrg( + operator, + init; + nsweeps=15, + maxdim=[10, 10, 20, 20, 40, 80, 100], + cutoff=(√(eps(real(elt)))), + noise=(√(eps(real(elt)))), + ) + rng = StableRNG(1234) + state = random_mps(rng, elt, s) + nexpansions = 10 + tau = elt(0.5) + for step in 1:nexpansions + # TODO: Use `fourthroot`/`∜` in Julia 1.10 and above. + state = expand( + state, operator; alg="global_krylov", krylovdim=3, cutoff=eps(real(elt))^(1//4) + ) + state = tdvp( + operator, + -4tau, + state; + nsteps=4, + cutoff=1e-5, + updater_kwargs=(; tol=1e-3, krylovdim=5), + ) + state = normalize(state) + end + @test scalartype(state) === elt + # TODO: Use `fourthroot`/`∜` in Julia 1.10 and above. + @test inner(state', operator, state) ≈ reference_energy rtol = 5 * eps(real(elt))^(1//4) + end +end +end diff --git a/test/base/test_solvers/test_exports.jl b/test/base/test_solvers/test_exports.jl new file mode 100644 index 0000000..ab457cd --- /dev/null +++ b/test/base/test_solvers/test_exports.jl @@ -0,0 +1,10 @@ +@eval module $(gensym()) +using ITensorTDVP: ITensorTDVP +using Test: @test, @testset +@testset "Test exports" begin + @test issetequal( + names(ITensorTDVP), + [:ITensorTDVP, :TimeDependentSum, :dmrg_x, :expand, :linsolve, :tdvp, :to_vec], + ) +end +end diff --git a/test/base/test_solvers/test_linsolve.jl b/test/base/test_solvers/test_linsolve.jl new file mode 100644 index 0000000..7079387 --- /dev/null +++ b/test/base/test_solvers/test_linsolve.jl @@ -0,0 +1,45 @@ +@eval module $(gensym()) +using ITensors: scalartype +using ITensors.ITensorMPS: MPO, OpSum, apply, random_mps, siteinds +using ITensorTDVP: ITensorTDVP, dmrg +using KrylovKit: linsolve +using LinearAlgebra: norm +using StableRNGs: StableRNG +using Test: @test, @test_throws, @testset +using Random: Random +@testset "linsolve (eltype=$elt, conserve_qns=$conserve_qns)" for elt in ( + Float32, Float64, Complex{Float32}, Complex{Float64} + ), + conserve_qns in [false, true] + + N = 6 + s = siteinds("S=1/2", N; conserve_qns) + os = OpSum() + for j in 1:(N - 1) + os += 0.5, "S+", j, "S-", j + 1 + os += 0.5, "S-", j, "S+", j + 1 + os += "Sz", j, "Sz", j + 1 + end + H = MPO(elt, os, s) + state = [isodd(n) ? "Up" : "Dn" for n in 1:N] + rng = StableRNG(1234) + x_c = random_mps(rng, elt, s, state; linkdims=2) + e, x_c = dmrg(H, x_c; nsweeps=10, cutoff=1e-6, maxdim=20, outputlevel=0) + @test scalartype(x_c) == elt + # Compute `b = H * x_c` + b = apply(H, x_c; cutoff=1e-8) + @test scalartype(b) == elt + # Starting guess + rng = StableRNG(1234) + x0 = x_c + elt(0.05) * random_mps(rng, elt, s, state; linkdims=2) + @test scalartype(x0) == elt + nsweeps = 10 + cutoff = 1e-5 + maxdim = 20 + updater_kwargs = (; tol=1e-4, maxiter=20, krylovdim=30, ishermitian=true) + @test_throws ErrorException linsolve(H, b, x0; cutoff, maxdim, updater_kwargs) + x = linsolve(H, b, x0; nsweeps, cutoff, maxdim, updater_kwargs) + @test scalartype(x) == elt + @test norm(x - x_c) < 1e-2 +end +end diff --git a/test/base/test_solvers/test_tdvp.jl b/test/base/test_solvers/test_tdvp.jl new file mode 100644 index 0000000..3e032d0 --- /dev/null +++ b/test/base/test_solvers/test_tdvp.jl @@ -0,0 +1,351 @@ +@eval module $(gensym()) +using ITensors: + ITensors, + AbstractObserver, + ITensor, + MPO, + MPS, + OpSum, + apply, + dag, + expect, + inner, + noprime, + op, + prime, + random_mps, + scalar, + siteinds +using ITensorTDVP: ITensorTDVP, tdvp +using KrylovKit: exponentiate +using LinearAlgebra: norm +using Observers: observer +using StableRNGs: StableRNG +using Test: @test, @test_throws, @testset +const elts = (Float32, Float64, Complex{Float32}, Complex{Float64}) +@testset "Basic TDVP (eltype=$elt)" for elt in elts + N = 10 + cutoff = eps(real(elt)) * 10^4 + s = siteinds("S=1/2", N) + os = OpSum() + for j in 1:(N - 1) + os += 0.5, "S+", j, "S-", j + 1 + os += 0.5, "S-", j, "S+", j + 1 + os += "Sz", j, "Sz", j + 1 + end + H = MPO(elt, os, s) + rng = StableRNG(1234) + ψ0 = random_mps(rng, elt, s; linkdims=10) + time_step = elt(0.1) * im + # Time evolve forward: + ψ1 = tdvp(H, -time_step, ψ0; cutoff, nsite=1) + @test ITensors.scalartype(ψ1) == complex(elt) + #Different backend updaters, default updater_backend = "exponentiate" + @test ψ1 ≈ tdvp(H, -time_step, ψ0; cutoff, nsite=1, updater_backend="applyexp") + @test norm(ψ1) ≈ 1 rtol = √eps(real(elt)) * 10 + ## Should lose fidelity: + #@test abs(inner(ψ0,ψ1)) < 0.9 + # Average energy should be conserved: + @test real(inner(ψ1', H, ψ1)) ≈ inner(ψ0', H, ψ0) rtol = √eps(real(elt)) * 10 + # Time evolve backwards: + ψ2 = tdvp(H, time_step, ψ1; cutoff) + @test ITensors.scalartype(ψ2) == complex(elt) + @test norm(ψ2) ≈ 1 rtol = √eps(real(elt)) * 10 + # Should rotate back to original state: + @test abs(inner(ψ0, ψ2)) > 0.99 +end + +@testset "TDVP: Sum of Hamiltonians (eltype=$elt)" for elt in elts + N = 10 + cutoff = 1e-10 + + s = siteinds("S=1/2", N) + + os1 = OpSum() + for j in 1:(N - 1) + os1 += 0.5, "S+", j, "S-", j + 1 + os1 += 0.5, "S-", j, "S+", j + 1 + end + os2 = OpSum() + for j in 1:(N - 1) + os2 += "Sz", j, "Sz", j + 1 + end + H1 = MPO(elt, os1, s) + H2 = MPO(elt, os2, s) + Hs = [H1, H2] + rng = StableRNG(1234) + ψ0 = random_mps(rng, elt, s; linkdims=10) + ψ1 = tdvp(Hs, -elt(0.1) * im, ψ0; cutoff, nsite=1) + @test ITensors.scalartype(ψ1) === complex(elt) + @test norm(ψ1) ≈ 1 rtol = √eps(real(elt)) + ## Should lose fidelity: + #@test abs(inner(ψ0,ψ1)) < 0.9 + # Average energy should be conserved: + @test real(sum(H -> inner(ψ1', H, ψ1), Hs)) ≈ sum(H -> inner(ψ0', H, ψ0), Hs) rtol = + 4 * √eps(real(elt)) + # Time evolve backwards: + ψ2 = tdvp(Hs, elt(0.1) * im, ψ1; cutoff) + @test ITensors.scalartype(ψ2) === complex(elt) + @test norm(ψ2) ≈ 1 rtol = √eps(real(elt)) + # Should rotate back to original state: + @test abs(inner(ψ0, ψ2)) > 0.99 +end +@testset "Custom updater in TDVP" begin + N = 10 + cutoff = 1e-12 + s = siteinds("S=1/2", N) + os = OpSum() + for j in 1:(N - 1) + os += 0.5, "S+", j, "S-", j + 1 + os += 0.5, "S-", j, "S+", j + 1 + os += "Sz", j, "Sz", j + 1 + end + H = MPO(os, s) + rng = StableRNG(1234) + ψ0 = random_mps(rng, s; linkdims=10) + function updater(PH, state0; internal_kwargs, kwargs...) + return exponentiate(PH, internal_kwargs.time_step, state0; kwargs...) + end + updater_kwargs = (; + ishermitian=true, tol=1e-12, krylovdim=30, maxiter=100, verbosity=0, eager=true + ) + t = -0.1im + ψ1 = tdvp(H, t, ψ0; updater, updater_kwargs, cutoff, nsite=1) + @test norm(ψ1) ≈ 1 + ## Should lose fidelity: + #@test abs(inner(ψ0,ψ1)) < 0.9 + # Average energy should be conserved: + @test real(inner(ψ1', H, ψ1)) ≈ inner(ψ0', H, ψ0) + # Time evolve backwards: + ψ2 = tdvp(H, +0.1im, ψ1; cutoff) + @test norm(ψ2) ≈ 1 + # Should rotate back to original state: + @test abs(inner(ψ0, ψ2)) > 0.99 +end +@testset "Accuracy Test" begin + N = 4 + tau = 0.1 + ttotal = 1.0 + cutoff = 1e-12 + s = siteinds("S=1/2", N; conserve_qns=false) + os = OpSum() + for j in 1:(N - 1) + os += 0.5, "S+", j, "S-", j + 1 + os += 0.5, "S-", j, "S+", j + 1 + os += "Sz", j, "Sz", j + 1 + end + H = MPO(os, s) + HM = prod(H) + Ut = exp(-im * tau * HM) + state = MPS(s, n -> isodd(n) ? "Up" : "Dn") + state2 = deepcopy(state) + statex = prod(state) + Sz_tdvp = Float64[] + Sz_tdvp2 = Float64[] + Sz_exact = Float64[] + c = div(N, 2) + Szc = op("Sz", s[c]) + Nsteps = Int(ttotal / tau) + for step in 1:Nsteps + statex = noprime(Ut * statex) + statex /= norm(statex) + + state = tdvp( + H, + -im * tau, + state; + cutoff, + normalize=false, + updater_kwargs=(; tol=1e-12, maxiter=500, krylovdim=25), + ) + push!(Sz_tdvp, real(expect(state, "Sz"; sites=c:c)[1])) + state2 = tdvp( + H, + -im * tau, + state2; + cutoff, + normalize=false, + updater_kwargs=(; tol=1e-12, maxiter=500, krylovdim=25), + ) + push!(Sz_tdvp2, real(expect(state2, "Sz"; sites=c:c)[1])) + push!(Sz_exact, real(scalar(dag(prime(statex, s[c])) * Szc * statex))) + F = abs(scalar(dag(statex) * prod(state))) + end + @test norm(Sz_tdvp - Sz_exact) < 1e-5 + @test norm(Sz_tdvp2 - Sz_exact) < 1e-5 +end +@testset "TEBD Comparison" begin + N = 10 + cutoff = 1e-12 + tau = 0.1 + ttotal = 1.0 + s = siteinds("S=1/2", N; conserve_qns=true) + os = OpSum() + for j in 1:(N - 1) + os += 0.5, "S+", j, "S-", j + 1 + os += 0.5, "S-", j, "S+", j + 1 + os += "Sz", j, "Sz", j + 1 + end + H = MPO(os, s) + gates = ITensor[] + for j in 1:(N - 1) + s1 = s[j] + s2 = s[j + 1] + hj = + op("Sz", s1) * op("Sz", s2) + + 1 / 2 * op("S+", s1) * op("S-", s2) + + 1 / 2 * op("S-", s1) * op("S+", s2) + Gj = exp(-1.0im * tau / 2 * hj) + push!(gates, Gj) + end + append!(gates, reverse(gates)) + state = MPS(s, n -> isodd(n) ? "Up" : "Dn") + phi = deepcopy(state) + c = div(N, 2) + # Evolve using TEBD + Nsteps = convert(Int, ceil(abs(ttotal / tau))) + Sz1 = zeros(Nsteps) + En1 = zeros(Nsteps) + Sz2 = zeros(Nsteps) + En2 = zeros(Nsteps) + for step in 1:Nsteps + state = apply(gates, state; cutoff) + nsite = (step <= 3 ? 2 : 1) + phi = tdvp( + H, -tau * im, phi; cutoff, nsite, normalize=true, updater_kwargs=(; krylovdim=15) + ) + Sz1[step] = expect(state, "Sz"; sites=c:c)[1] + Sz2[step] = expect(phi, "Sz"; sites=c:c)[1] + En1[step] = real(inner(state', H, state)) + En2[step] = real(inner(phi', H, phi)) + end + # Evolve using TDVP + function measure_sz(; state, bond, half_sweep) + if bond == 1 && half_sweep == 2 + return expect(state, "Sz"; sites=c) + end + return nothing + end + function measure_en(; state, bond, half_sweep) + if bond == 1 && half_sweep == 2 + return real(inner(state', H, state)) + end + return nothing + end + obs = observer("Sz" => measure_sz, "En" => measure_en) + + phi = MPS(s, n -> isodd(n) ? "Up" : "Dn") + phi = tdvp( + H, -im * ttotal, phi; time_step=-im * tau, cutoff, normalize=false, (observer!)=obs + ) + Sz2 = obs.Sz + En2 = obs.En + @test norm(Sz1 - Sz2) < 1e-3 + @test norm(En1 - En2) < 1e-3 +end +@testset "Imaginary Time Evolution" for reverse_step in [true, false] + N = 10 + cutoff = 1e-12 + tau = 1.0 + ttotal = 50.0 + s = siteinds("S=1/2", N) + os = OpSum() + for j in 1:(N - 1) + os += 0.5, "S+", j, "S-", j + 1 + os += 0.5, "S-", j, "S+", j + 1 + os += "Sz", j, "Sz", j + 1 + end + H = MPO(os, s) + rng = StableRNG(1234) + state = random_mps(rng, s; linkdims=2) + state2 = deepcopy(state) + trange = 0.0:tau:ttotal + for (step, t) in enumerate(trange) + nsite = (step <= 10 ? 2 : 1) + state = tdvp( + H, + -tau, + state; + cutoff, + nsite, + reverse_step, + normalize=true, + updater_kwargs=(; krylovdim=15), + ) + state2 = tdvp( + H, + -tau, + state2; + cutoff, + nsite, + reverse_step, + normalize=true, + updater_kwargs=(; krylovdim=15), + ) + end + @test state ≈ state2 rtol = 1e-6 + en1 = inner(state', H, state) + en2 = inner(state2', H, state2) + @test en1 < -4.25 + @test en1 ≈ en2 +end +@testset "Observers" begin + N = 10 + cutoff = 1e-12 + tau = 0.1 + ttotal = 1.0 + s = siteinds("S=1/2", N; conserve_qns=true) + os = OpSum() + for j in 1:(N - 1) + os += 0.5, "S+", j, "S-", j + 1 + os += 0.5, "S-", j, "S+", j + 1 + os += "Sz", j, "Sz", j + 1 + end + H = MPO(os, s) + c = div(N, 2) + # Using Observers.jl + function measure_sz(; state, bond, half_sweep) + if bond == 1 && half_sweep == 2 + return expect(state, "Sz"; sites=c) + end + return nothing + end + function measure_en(; state, bond, half_sweep) + if bond == 1 && half_sweep == 2 + return real(inner(state', H, state)) + end + return nothing + end + function identity_info(; info) + return info + end + obs = observer("Sz" => measure_sz, "En" => measure_en, "info" => identity_info) + step_measure_sz(; state) = expect(state, "Sz"; sites=c) + step_measure_en(; state) = real(inner(state', H, state)) + step_obs = observer("Sz" => step_measure_sz, "En" => step_measure_en) + state = MPS(s, n -> isodd(n) ? "Up" : "Dn") + tdvp( + H, + -im * ttotal, + state; + time_step=-im * tau, + cutoff, + normalize=false, + (observer!)=obs, + (step_observer!)=step_obs, + ) + Sz = filter(!isnothing, obs.Sz) + En = filter(!isnothing, obs.En) + infos = obs.info + Sz_step = step_obs.Sz + En_step = step_obs.En + @test length(Sz) == 10 + @test length(En) == 10 + @test length(Sz_step) == 10 + @test length(En_step) == 10 + @test Sz ≈ Sz_step + @test En ≈ En_step + @test all(x -> x.info.converged == 1, infos) + @test length(values(infos)) == 180 +end +end diff --git a/test/base/test_solvers/test_tdvp_time_dependent.jl b/test/base/test_solvers/test_tdvp_time_dependent.jl new file mode 100644 index 0000000..26b3290 --- /dev/null +++ b/test/base/test_solvers/test_tdvp_time_dependent.jl @@ -0,0 +1,109 @@ +@eval module $(gensym()) +using ITensors: ITensors, Index, QN, contract, scalartype +using ITensors.ITensorMPS: MPO, MPS, ProjMPO, ProjMPOSum, random_mps, position!, siteinds +using ITensorTDVP: ITensorTDVP, TimeDependentSum, tdvp +using LinearAlgebra: norm +using StableRNGs: StableRNG +using Test: @test, @test_skip, @testset +include(joinpath(pkgdir(ITensorTDVP), "examples", "03_models.jl")) +include(joinpath(pkgdir(ITensorTDVP), "examples", "03_updaters.jl")) +@testset "TDVP with ODE local updater" begin + @testset "TimeDependentSum (eltype=$elt)" for elt in ( + Float32, Float64, Complex{Float32}, Complex{Float64} + ), + conserve_qns in [false, true] + + n = 4 + s = siteinds("S=1/2", 4; conserve_qns) + H = MPO(elt, s, "I") + H⃗ = (H, H) + region = 2:3 + rng = StableRNG(1234) + ψ = random_mps(rng, elt, s, j -> isodd(j) ? "↑" : "↓"; linkdims=2) + H⃗ᵣ = ProjMPO.(H⃗) + map(Hᵣ -> position!(Hᵣ, ψ, first(region)), H⃗ᵣ) + ∑Hᵣ = ProjMPOSum(collect(H⃗)) + position!(∑Hᵣ, ψ, first(region)) + f⃗ₜ = (t -> sin(elt(0.1) * t), t -> cos(elt(0.2) * t)) + α = elt(0.5) + ∑Hₜ = α * TimeDependentSum(f⃗ₜ, ITensors.terms(∑Hᵣ)) + t₀ = elt(0.5) + ∑Hₜ₀ = ∑Hₜ(t₀) + ψᵣ = reduce(*, map(v -> ψ[v], region)) + Hψ = ∑Hₜ₀(ψᵣ) + @test eltype(Hψ) == elt + @test Hψ ≈ sum(i -> α * f⃗ₜ[i](t₀) * H⃗ᵣ[i](ψᵣ), eachindex(H⃗)) + end + @testset "Time dependent Hamiltonian (eltype=$elt, conserve_qns=$conserve_qns)" for elt in + ( + Float32, Float64, Complex{Float32}, Complex{Float64} + ), + conserve_qns in [false, true] + + n = 4 + J₁ = elt(1) + J₂ = elt(0.1) + ω₁ = real(elt)(0.1) + ω₂ = real(elt)(0.2) + ω⃗ = (ω₁, ω₂) + f⃗ = map(ω -> (t -> cos(ω * t)), ω⃗) + time_step = real(elt)(0.1) + time_stop = real(elt)(1) + nsite = 2 + maxdim = 100 + cutoff = √(eps(real(elt))) + tol = √eps(real(elt)) + s = siteinds("S=1/2", n) + ℋ₁₀ = heisenberg(n; J=J₁, J2=zero(elt)) + ℋ₂₀ = heisenberg(n; J=zero(elt), J2=J₂) + ℋ⃗₀ = (ℋ₁₀, ℋ₂₀) + H⃗₀ = map(ℋ₀ -> MPO(elt, ℋ₀, s), ℋ⃗₀) + ψ₀ = complex.(MPS(elt, s, j -> isodd(j) ? "↑" : "↓")) + ψₜ_ode = tdvp( + -im * TimeDependentSum(f⃗, H⃗₀), + time_stop, + ψ₀; + updater=ode_updater, + updater_kwargs=(; reltol=tol, abstol=tol), + time_step, + maxdim, + cutoff, + nsite, + ) + ψₜ_krylov = tdvp( + -im * TimeDependentSum(f⃗, H⃗₀), + time_stop, + ψ₀; + updater=krylov_updater, + updater_kwargs=(; tol, eager=true), + time_step, + maxdim, + cutoff, + nsite, + ) + ψₜ_full, _ = ode_updater( + -im * TimeDependentSum(f⃗, contract.(H⃗₀)), + contract(ψ₀); + internal_kwargs=(; time_step=time_stop), + reltol=tol, + abstol=tol, + ) + + @test scalartype(ψ₀) == complex(elt) + @test scalartype(ψₜ_ode) == complex(elt) + @test scalartype(ψₜ_krylov) == complex(elt) + @test scalartype(ψₜ_full) == complex(elt) + @test norm(ψ₀) ≈ 1 + @test norm(ψₜ_ode) ≈ 1 + @test norm(ψₜ_krylov) ≈ 1 rtol = √(eps(real(elt))) + @test norm(ψₜ_full) ≈ 1 + + ode_err = norm(contract(ψₜ_ode) - ψₜ_full) + krylov_err = norm(contract(ψₜ_krylov) - ψₜ_full) + + @test krylov_err > ode_err + @test ode_err < √(eps(real(elt))) * 10^4 + @test krylov_err < √(eps(real(elt))) * 10^5 + end +end +end diff --git a/test/backup/utils/TestITensorMPSExportedNames.jl b/test/base/utils/TestITensorMPSExportedNames.jl similarity index 100% rename from test/backup/utils/TestITensorMPSExportedNames.jl rename to test/base/utils/TestITensorMPSExportedNames.jl diff --git a/test/ext/ITensorMPSChainRulesCoreExt/test_chainrules.jl b/test/ext/ITensorMPSChainRulesCoreExt/test_chainrules.jl index 9b5e493..c704920 100644 --- a/test/ext/ITensorMPSChainRulesCoreExt/test_chainrules.jl +++ b/test/ext/ITensorMPSChainRulesCoreExt/test_chainrules.jl @@ -1,3 +1,4 @@ +using ITensorMPS using ITensors using Random using Test