diff --git a/Project.toml b/Project.toml index 4a658ef..e2ce40a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "TeNLib" uuid = "1ce00845-9e1b-4798-bd63-033d15ad4ca9" authors = ["Titas Chanda"] -version = "0.0.1-DEV" +version = "0.1.0" [deps] DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" diff --git a/docs/make.jl b/docs/make.jl index 2f0e244..1892ae4 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -34,7 +34,12 @@ settings = Dict( "Tree Tensor Networks" => "ttn/ttn.md", "The TTN object" => "ttn/ttn_struct.md", "Generating TTN" => "ttn/ttn_gen.md", - "StateEnvsTTN" => "ttn/state_envs_ttn.md" + "StateEnvsTTN" => "ttn/state_envs_ttn.md", + "Perform local updates" => "ttn/update_site_ttn.md", + "Sweeping through the TTN" => "ttn/sweep_ttn.md", + "Optimizing TTN" => "ttn/optimize_ttn.md", + "Example: Optimizing TTN" => "ttn/example_optimize.md", + "Measuring the TTN" => "ttn/measure_ttn.md" ] ], :format => Documenter.HTML( diff --git a/docs/src/base/couplingmodel.md b/docs/src/base/couplingmodel.md index d16820c..122f31e 100644 --- a/docs/src/base/couplingmodel.md +++ b/docs/src/base/couplingmodel.md @@ -1,7 +1,7 @@ # `CouplingModel` TeNLib.jl degines as struct, called the `CouplingModel`, to store the Hamiltonian terms. -In case of MPS based algorithms, `CouplingModel` can be replaced by `MPO` without modifying +In case of MPS based algorithms, `CouplingModel` can replace `MPO` without modifying rest of the code. For Tree Tensor Network (TTN) codes, only `CouplingModel` can be used. Different elements of `CouplingModel` are contracted in parallel. diff --git a/docs/src/index.md b/docs/src/index.md index 3a3ca80..8916a29 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -14,7 +14,7 @@ The documentation for TeNLib.jl can be found [**here**](https://titaschanda.gith TeNLib.jl features widely-used Tensor Network (TN) codes, designed with a **multi-layered abstraction** to cater to diverse user needs. The library provides users with varying levels of control over their computations. -Presently, TeNLib.jl presents an array of functionalities for: +Currently, TeNLib.jl presents an array of functionalities for: * *(a)* **Finite-size Matrix-Product States (MPS)**: Different variants of Density Matrix Renormalization Group (DMRG) and Time Dependent Variational Principle (TDVP) (including subspace expansion) methods. * *(b)* **Tree Tensor Network (TTN)**: Variational search for the ground state and first few excited states. @@ -33,6 +33,17 @@ pkg> add ITensors pkg> add https://github.com/titaschanda/TeNLib.jl ``` +## Found an Issue or Bug? + +> "Beware of bugs in the above code; I have only proved it correct, not tried it." +> -- Donald Knuth + + +If you find bugs or a mistakes of any kind, please let us know by adding an issue to the +[GitHub issue tracker](https://github.com/titaschanda/TeNLib.jl/issues). +You are also welcome to submit a [pull request](https://github.com/titaschanda/TeNLib.jl/pulls). + + ## Future functionality? Here is a list for future additions in the decreasing order of priority. Any help / suggestion is welcome. @@ -41,6 +52,13 @@ Here is a list for future additions in the decreasing order of priority. Any hel * **Projected Entangled Pair States (PEPS)** for 2D problems. * Real-time evolution method using PEPS and TTN. +Also, please feel free to ask about a new feature by adding a new request to the +[GitHub issue tracker](https://github.com/titaschanda/TeNLib.jl/issues) labelled +`feature request`. Note that submitting a pull request, providing the needed changes to +introduced your requested feature, will speed up the process. + + + ## Example: A simple DMRG code The following code is for a simple DMRG run at **the highest level of abstraction without any additional control**. @@ -145,7 +163,4 @@ let en, psi = optimize(psi0, H, params, sweeppath) end -``` -!!! info - "`OpStrings` and `CouplingModel` can be also used for MPS based codes without modifying other parts of the code." - \ No newline at end of file +``` \ No newline at end of file diff --git a/docs/src/mps/example_dmrg.md b/docs/src/mps/example_dmrg.md index 042a202..56ddd60 100644 --- a/docs/src/mps/example_dmrg.md +++ b/docs/src/mps/example_dmrg.md @@ -59,7 +59,7 @@ psi = getpsi(sysenv) # psi = sysenv.psi ``` Most often, it is better to do a single-site DMRG (without any noise) after standard two-site update for better -convergence. Such lowe-level function using `StateEnvs` is useful for that. +convergence. Such lower-level function using `StateEnvs` is useful for that. ``` sysenv = StateEnvs(psi0, H) @@ -83,8 +83,9 @@ built from a single MPO. ``` krylov_extend!(sysenv; extension_applyH_maxdim = 40) ``` -**Note**: Global Subspace Expansion can result into huge MPS bond dimension. That is why -the named input parameters of [`krylov_extend!`](@ref krylov_extend!(sysenv::StateEnvs{ProjMPO}; kwargs...)) should be chosen carefully. +!!! warning + Global Subspace Expansion may result into huge MPS bond dimension. That is why + the named input parameters of [`krylov_extend!`](@ref krylov_extend!(sysenv::StateEnvs{ProjMPO}; kwargs...)) should be chosen carefully. ## Excited state DMRG diff --git a/docs/src/mps/sweep.md b/docs/src/mps/sweep.md index ae75952..17cbc46 100644 --- a/docs/src/mps/sweep.md +++ b/docs/src/mps/sweep.md @@ -31,14 +31,15 @@ dynamic_fullsweep!(sysenv::StateEnvs, solver, swdata::SweepData; kwargs...) ## Global Subspace Expansion -Following [Phys. Rev. B **102**, 094315 (2020)](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.102.094315), a Global Subspace Expansion can be performed using Krylov subspace if the -environments are created by a single `MPO`. +Following [Phys. Rev. B **102**, 094315 (2020)](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.102.094315), a Global Subspace Expansion can be performed using Krylov subspace if the environments are created by a single `MPO`. -Apart from TDVP, Global Subspace Expansion is also very useful for DMRG to get rid of nasty -local minimas. ```@docs krylov_extend!(psi::MPS, H::MPO; kwargs...) krylov_extend!(sysenv::StateEnvs{ProjMPO}; kwargs...) ``` +!!! info + Apart from TDVP, Global Subspace Expansion is also very useful for DMRG to get rid of nasty + local minimas. + diff --git a/docs/src/ttn/example_optimize.md b/docs/src/ttn/example_optimize.md new file mode 100644 index 0000000..45d8e83 --- /dev/null +++ b/docs/src/ttn/example_optimize.md @@ -0,0 +1,88 @@ +# Example: Optimizing TTN + +## Ground state search + +A straight-forward TTN optimization can be performed as follows. + +``` +using ITensors +using TeNLib + +let + N = 32 + sites = siteinds("S=1/2",N) + os = OpStrings() + + for j=1:N-1 + os += 1, "Sz" => j, "Sz" => j+1 + os += 0.5, "S+" => j, "S-" => j+1 + os += 0.5, "S-" => j, "S+" => j+1 + end + + H = CouplingModel(os,sites) + + psi0 = default_randomTTN(sites, 12, QN("Sz", 0)) + sweeppath = default_sweeppath(psi0) + + params = OptimizeParamsTTN(; nsweeps = [10, 10], maxdim = [20, 50], + cutoff = 1e-14, noise = 1e-3, noisedecay = 2, + disable_noise_after = 5) + + en, psi = optimize(psi0, H, params, sweeppath) +end +``` + +!!! warning + For TTN optimzaion, direct `MPO` input is not supported. One can, however, prepare a + `CouplingModel` from `MPO` using + [`CouplingModel(mpos::MPO...)`](@ref CouplingModel(mpos::MPO...)). + + +Instead of using such higher-level code, one can also use lower-level functions for a better +control. +``` +sysenv = StateEnvsTTN(psi0, H) + +swdata = optimize!(sysenv, params, sweeppath) + +# Get energy from `Sweepdata` +energy = swdata.energy[end] + +# take a shallow copy of the TTN +# if the `StateEnvsTTN` will be updated later again +psi = getpsi(sysenv) + +# Alternatively, take the psi from `StateEnvsTTN` itself. +# NOTE: This can crash the simulation, if the TTN is modified (e.g., in measurements) +# and `StateEnvsTTN` is going to be updated later on. +# psi = sysenv.psi +``` + +## Excited state search + +Excited state search is also straightforword. + +``` +# Given a ground state `psi_gr`, initial TTN `psi0`, +# and a Hamiltonian `H` + +en, psi = optimize(psi0, H, [psi_gr], params, sweeppath; weight = 10.0) +``` + +Similarly, using `StateEnvsTTN`: +``` +sysenv_ex = StateEnvsTTN(psi0, H, [psi_gr]; weight = 10.0) +swdata_ex = optimize!(sysenv_ex, params, sweeppath) + +# Get energy from `Sweepdata` +energy1 = swdata_ex.energy[end] + +# take a shallow copy of the TTN +# if the `StateEnvsTTN` will be updated later again +psi1 = getpsi(sysenv_ex) + +# Alternatively, take the psi from `StateEnvsTTN` itself. +# NOTE: This can crash the simulation, if the TTN is modified (e.g., in measurements) +# and `StateEnvsTTN` is going to be updated later. +# psi1 = sysenv_ex.psi +``` \ No newline at end of file diff --git a/docs/src/ttn/measure_ttn.md b/docs/src/ttn/measure_ttn.md new file mode 100644 index 0000000..83dddae --- /dev/null +++ b/docs/src/ttn/measure_ttn.md @@ -0,0 +1,9 @@ +# Measuring the TTN + +```@docs +measure(::Type{T}, psi::TTN, opten::ITensor) where T <: Union{ComplexF64, Float64} +measure(::Type{T}, psi::TTN, opstr::String, pos::Int) where T <: Union{ComplexF64, Float64} +measure(::Type{T}, psi::TTN, opstr::String; kwargs...) where T <: Union{ComplexF64, Float64} +measure(::Type{T}, psi::TTN, optens::Vector{ITensor}) where T <: Union{ComplexF64, Float64} +measure(::Type{T}, psi::TTN, oppairs::Vector{Pair{String, Int}}; isfermions::Bool = true) where T <: Union{ComplexF64, Float64} +``` \ No newline at end of file diff --git a/docs/src/ttn/optimize_ttn.md b/docs/src/ttn/optimize_ttn.md new file mode 100644 index 0000000..a8e5a84 --- /dev/null +++ b/docs/src/ttn/optimize_ttn.md @@ -0,0 +1,29 @@ +# Optimizing TTN + +Functions to optimize TTN. + +## `OptimizeParamsTTN` + +TeNLib.jl defines a struct, called `OptimizeParamsTTN`, to control TTN optimizations. + +```@docs +OptimizeParamsTTN +OptimizeParamsTTN(;maxdim::Vector{Int}, nsweeps::Vector{Int}, cutoff::Union{Vector{Float64}, Float64} = _Float64_Threshold, noise::Union{Vector{Float64}, Float64, Int} = 0.0, noisedecay::Union{Vector{Float64}, Float64, Int} = 1.0, disable_noise_after::Union{Vector{Int}, Int} = typemax(Int)) +Base.copy(params::OptimizeParamsTTN) +``` + +## A lower level optimization function + +Following function modifies `StateEnvsTTN` in-place. Skip this function if you want to +avoid lower-level abstraction. + +```@docs +optimize!(sysenv::StateEnvsTTN, params::OptimizeParamsTTN,sweeppath::Vector{Int2}; kwargs...) +``` + +## Higher level optimzation functions + +```@docs +optimize(psi0::TTN, H::CouplingModel, params::OptimizeParamsTTN, sweeppath::Vector{Int2}; kwargs...) +``` + diff --git a/docs/src/ttn/sweep_ttn.md b/docs/src/ttn/sweep_ttn.md new file mode 100644 index 0000000..fff52e2 --- /dev/null +++ b/docs/src/ttn/sweep_ttn.md @@ -0,0 +1,26 @@ +# Sweeping through the TTN + +At a lower level of abstraction, TeNLib.jl allows to control each fullsweep + manually to update `StateEnvsTTN`. + +Skip this part if you want to avoid lower-level abstraction. + +## `SweepDataTTN` + +TeNLib.jl defines a struct, called `SweepDataTTN`, to store essential data after each fullsweep. + +```@docs +SweepDataTTN +Base.copy(swdata::SweepDataTTN) +``` +## `sweeppath` + +```@docs +default_sweeppath +``` + +## Perform a fullsweep + +```@docs +fullsweep!(sysenv::StateEnvsTTN, sweeppath::Vector{Int2}, solver, swdata::SweepDataTTN; kwargs...) +``` diff --git a/docs/src/ttn/ttn.md b/docs/src/ttn/ttn.md index ad644e4..79dfd5a 100644 --- a/docs/src/ttn/ttn.md +++ b/docs/src/ttn/ttn.md @@ -21,11 +21,18 @@ described in [Phys. Rev. B **90**, 125154 (2014)](https://journals.aps.org/prb/a --- Each tensor in the TTN are indexed by a pair (`Tuple`) of `Int`s = `(ll, nn)`. -``` -const Int2 = Tuple{Int, Int} +```@docs +Int2 ``` In the default scenario, `ll` denotes the layer index and `nn` denotes the tensor index at each layer (see the image above). The counting for `ll`, in the default case, starts at the bottom (towards to top level), while `nn` counts from left. The structure of the network is defined by a -[`Graph{Int2}`](@ref "The Graph object") object. +[`Graph{Int2}`](@ref "The Graph object") object. The link / bond between two neighboring nodes is +denoted by a `LinkTypeTTN` object, an **unordered** pair of `Int2`s. +```@docs +LinkTypeTTN +``` +!!! info + The order of the nodes inside `LinkTypeTTN` is irrelevant, i.e., + `LinkTypeTTN(node1, node2) == LinkTypeTTN(node2, node1)`. --- diff --git a/docs/src/ttn/update_site_ttn.md b/docs/src/ttn/update_site_ttn.md new file mode 100644 index 0000000..e66b7fc --- /dev/null +++ b/docs/src/ttn/update_site_ttn.md @@ -0,0 +1,16 @@ +# Perform local updates + +At the lowest-level of abstraction, TeNLib.jl allows for updating the `StateEnvsTTN` for each sites manually. + +Skip this part if you want to avoid lower-level abstraction. + +```@docs +update_position!(sysenv::StateEnvsTTN, solver, node::Int2; time_step::Union{Float64, ComplexF64, Nothing}, normalize::Bool, maxdim::Int, mindim::Int, cutoff::Float64, svd_alg::String, kwargs...) +``` + +TeNLib.jl implements the subspace expansion method described in +[SciPost Phys. Lect. Notes 8 (2019)] (https://scipost.org/10.21468/SciPostPhysLectNotes.8) to increase the bond dimension between two neighboring nodes. + +```@docs +subspace_expand!(psi::TTN, node::Int2, nextnode::Int2, max_expand_dim::Int, noise::Float64) +``` \ No newline at end of file diff --git a/src/export.jl b/src/export.jl index 99d8a29..637aa56 100644 --- a/src/export.jl +++ b/src/export.jl @@ -3,9 +3,7 @@ export # base/global_variables.jl - _Float64_threshold, Float64_threshold, - _using_threaded_loop, using_threaded_loop, enable_threaded_loop, disable_threaded_loop, @@ -15,13 +13,9 @@ export Vector2, IDType, IDTensors, - eltype, # base/helper_internal_funcs.jl gen_rand_id, - _divide_by_chunksize, - _add_oplinks, - _directsum, combineinds, indexintersection, @@ -43,7 +37,6 @@ export nodes_from_bfs, shortest_path, nextnode_in_path, - _has_cycle_dfs, has_cycle, find_sum_central_node, find_eccentric_central_node, @@ -63,23 +56,16 @@ export mergeterms, # base/opstringsmpo.jl - _chunksum_mpoterms, # base/couplingmodel.jl CouplingModel, - _initCouplingModel, # mps/projmps2.jl ProjMPS2, set_nsite!, - _makeL!, - makeL!, - _makeR!, - makeR!, position!, contract, - proj_mps, product, # mps/projmpo_mps2.jl @@ -108,15 +94,7 @@ export ProjCouplingModel, set_nsite!, nsite, - site_range, - lproj, - rproj, - _makeL!, - makeL!, - _makeR!, - makeR!, position, - _contract, product, # mps/projcouplingmodel_mps.jl @@ -137,9 +115,6 @@ export product, # base/update_site.jl - halfsweep_done, - _update_two_site!, - _update_one_site!, update_position!, # mps/sweep.jl @@ -147,7 +122,6 @@ export fullsweep!, dynamic_fullsweep!, krylov_extend!, - _krylov_addbasis, # mps/dmrg.jl DMRGParams, @@ -196,21 +170,14 @@ export isometrize!, # ttn/helper_internal_funcs.jl - _minimum_power2_greater_than, - _get_links, # ttn/ttn_generators.jl - _distribute_site_positions, default_graph_sitenodes, - _ttn_ind_reducedim!, - _ttn_ind_cleanup_one!, - _ttn_ind_cleanup!, randomTTN, default_randomTTN, # ttn/linktensors.jl LinkTensorsTTN, - _collect_link_tensors, move_linktensors_to_next!, move_linktensors!, product, @@ -225,8 +192,7 @@ export LinkProjTTN, # ttn/measure_ttn.jl - expectC, - expectR, + measure, # ttn/environment.jl EnvCouplingModelTTN, diff --git a/src/mps/dmrg.jl b/src/mps/dmrg.jl index 171aa44..34aa7b1 100644 --- a/src/mps/dmrg.jl +++ b/src/mps/dmrg.jl @@ -134,7 +134,7 @@ to trigger this early stopping. #### Named arguments for `solver` and their default values: See the documentation of KrylovKit.jl. - - `ishermitian::Bool = true` + - `ishermitian::Bool = true`. - `solver_tol::Float64 = 1E-14`. - `solver_krylovdim::Int = 5`. - `solver_maxiter::Int = 2`. @@ -280,7 +280,7 @@ to trigger this early stopping. #### Named arguments for `eig_solver` and their default values: See the documentation of KrylovKit.jl. - - `ishermitian::Bool = true` + - `ishermitian::Bool = true`. - `solver_tol::Float64 = 1E-14`. - `solver_krylovdim::Int = 5`. - `solver_maxiter::Int = 2`. diff --git a/src/mps/sweep.jl b/src/mps/sweep.jl index 1df0d89..cb625ce 100644 --- a/src/mps/sweep.jl +++ b/src/mps/sweep.jl @@ -21,8 +21,7 @@ Holds historical data after each (full)sweep. Requires for convergence check etc previous halfsweep. #### Default constructor: - - `SweepData()`: Initialize an empty `SweepData` object. Required to call following functions -at the first time. + - `SweepData()`: Initialize an empty `SweepData` object. """ mutable struct SweepData sweepcount::Int @@ -42,7 +41,7 @@ SweepData() = SweepData(0, Int[], Float64[], Float64[], Float64[], Float64[]) Shallow copy of `SweepData`. """ Base.copy(swdata::SweepData) = SweepData(swdata.sweepcount, - Base.copy(swdata.sweepcount), + Base.copy(swdata.maxchi), Base.copy(swdata.energy), Base.copy(swdata.entropy), Base.copy(swdata.mastruncerr), @@ -78,7 +77,7 @@ Perform a fullsweep (left-to-right and right-to-left) by `solver`. #### Named arguments for `solver` and their default values: See the documentation of KrylovKit.jl. - - `ishermitian::Bool = true` + - `ishermitian::Bool = true`. - `solver_tol::Float64 = 1E-14` if `eig_solver`, `1E-12` if `exp_solver`. - `solver_krylovdim::Int = 5` if `eig_solver`, `30` if `exp_solver`. - `solver_maxiter::Int = 2` if `eig_solver`, `100` if `exp_solver`. @@ -226,7 +225,7 @@ that bond in the subsequent halfsweep, otherwise performs two-site update. #### Named arguments for `solver` and their default values: See the documentation of KrylovKit.jl. - - `ishermitian::Bool = true` + - `ishermitian::Bool = true`. - `solver_tol::Float64 = 1E-14` if `eig_solver`, `1E-12` if `exp_solver`. - `solver_krylovdim::Int = 3` if `eig_solver`, `30` if `exp_solver`. - `solver_maxiter::Int = 1` if `eig_solver`, `100` if `exp_solver`. diff --git a/src/mps/tdvp.jl b/src/mps/tdvp.jl index 75a41f5..b4bd1d3 100644 --- a/src/mps/tdvp.jl +++ b/src/mps/tdvp.jl @@ -227,7 +227,7 @@ Therefore, for real-time dynamics with step `dt`, `time_step` should be `-im * d #### Named arguments for `solver` and their default values: See the documentation of KrylovKit.jl. - - `ishermitian::Bool = true` + - `ishermitian::Bool = true`. - `solver_tol::Float64 = 1E-12`. - `solver_krylovdim::Int = 30`. - `solver_maxiter::Int = 100`. diff --git a/src/mps/update_site.jl b/src/mps/update_site.jl index 4c43da4..a955f32 100644 --- a/src/mps/update_site.jl +++ b/src/mps/update_site.jl @@ -176,7 +176,7 @@ end """ update_position!(sysenv::StateEnvs, solver, pos::Int, nsite::Int, ortho::String; kwargs...) -Update StateEnvs at position `pos` by `solver`. +Updates StateEnvs at position `pos` by `solver`. #### Arguments: - `sysenv::StateEnvs` @@ -197,7 +197,7 @@ Update StateEnvs at position `pos` by `solver`. #### Named arguments for `solver` and their default values: See the documentation of KrylovKit.jl. - - `ishermitian::Bool = true` + - `ishermitian::Bool = true`. - `solver_tol::Float64 = 1E-14` if `eig_solver`, `1E-12` if `exp_solver`. - `solver_krylovdim::Int = 5` if `eig_solver`, `30` if `exp_solver`. - `solver_maxiter::Int = 2` if `eig_solver`, `100` if `exp_solver`. @@ -208,7 +208,7 @@ See the documentation of KrylovKit.jl. #### Return values: - `::Float64`: Energy. - `::Float64`: Truncation Error. - - `Vector{Float64}`: SVD spectrum. + - `::Vector{Float64}`: SVD spectrum. """ function update_position!(sysenv::StateEnvs, solver, pos::Int, nsite::Int, ortho::String; diff --git a/src/ttn/linktensors.jl b/src/ttn/linktensors.jl index d56a4a1..5add4da 100644 --- a/src/ttn/linktensors.jl +++ b/src/ttn/linktensors.jl @@ -301,7 +301,7 @@ function LinkTensorsTTN(psi::TTN, optens::Vector{ITensor})::LinkTensorsTTN positions = sort(collect(keys(opdict))) for ii = 2 : length(positions) - dummyindex = hasqns(sites) ? + dummyindex = !removeqn ? Index(QN() => 1; tags = "OpLink") : Index(1; tags = "OpLink") opdict[positions[ii-1]] *= onehot(dummyindex => 1) diff --git a/src/ttn/measure_ttn.jl b/src/ttn/measure_ttn.jl index a706160..97d99ab 100644 --- a/src/ttn/measure_ttn.jl +++ b/src/ttn/measure_ttn.jl @@ -1,7 +1,16 @@ ################################################################################# -function expectC(psi::TTN, opten::ITensor)::ComplexF64 +""" + function measure(::Type{ElT} = ComplexF64, psi::TTN, opten::ITensor) + +Returns (complex / real) local expectation value (`::ComplexF64` / `::Float64`) +of a given TTN `psi::TTN`. The operator `opten::ITensor` must be single-site operator. + +For `ElT = Float64`, if the expectation value is complex, raises a warning! +""" +function measure(::Type{T}, psi::TTN, + opten::ITensor)::T where T <: Union{ComplexF64, Float64} pos = findsite(psi, opten) pass = hasind(opten, siteind(psi, pos)) && hasind(opten, siteind(psi, pos)') && @@ -17,47 +26,75 @@ function expectC(psi::TTN, opten::ITensor)::ComplexF64 ket = conditional_removeqns(psi[node]) ind = commonind(ket, opten) bra = dag(prime(ket, ind)) - val = scalar(bra * opten * ket) - return complex(val) -end + val = complex(scalar(bra * opten * ket)) -################################################################################# - -function expectR(psi::TTN, opten::ITensor)::Float64 - val = expectC(psi, opten) - if abs(imag(val)) > 100 * Float64_threshold() - @warn "`expectR()`: Got complex number with |imaginary part| > $(100 * Float64_threshold()), use `expectC instead !!" + if T <: Complex + return val + else + if abs(imag(val)) > 100 * Float64_threshold() + @warn "`measure(::Float64)`: Got complex number with |imaginary part| > $(100 * Float64_threshold()) !!" + end + return real(val) end - return real(val) + + return complex(val) end +measure(psi::TTN, opten::ITensor) = measure(ComplexF64, psi, opten) + ################################################################################# -expectC(psi::TTN, opstr::String, pos::Int) = expectC(psi, op(opstr, siteind(psi, pos))) +""" + function measure(::Type{ElT} = ComplexF64, psi::TTN, opstr::String, pos::Int) -################################################################################# +Returns (complex / real) local expectation value (`::ComplexF64` / `::Float64`) of a +given TTN `psi::TTN` for a given operator name (`opstr::String`) and a site index (`pos::Int`) + +For `ElT = Float64`, if the expectation value is complex, raises a warning! +""" +measure(::Type{T}, psi::TTN, opstr::String, + pos::Int) where T <: Union{ComplexF64, Float64} = + measure(T, psi, op(opstr, siteind(psi, pos))) -expectR(psi::TTN, opstr::String, pos::Int) = expectR(psi, op(opstr, siteind(psi, pos))) +measure(psi::TTN, opstr::String, pos::Int) = + measure(ComplexF, psi, opstr, Int) ################################################################################# -function expectC(psi::TTN, opstr::String; kwargs...) - N = numsites(psi) - sites = get(kwargs, :sites, 1:N) - return map(pos -> expectC(psi, opstr, pos), sites) -end +""" + function measure(::Type{ElT} = ComplexF64, psi::TTN, opstr::String; kwargs...) -################################################################################# +Returns (complex / real) local expectation values (`::Vector{ComplexF64}` / `::Vector{Float64}`) +of a given TTN `psi::TTN` for a given operator name (`opstr::String`) at all the sites. + +Optionally, for specific sites, keyword argument `sites` can be specified, e.g., +`sites = [1, 2, 3]`. -function expectR(psi::TTN, opstr::String; kwargs...) +For `ElT = Float64`, if the expectation value is complex, raises a warning! +""" +function measure(::Type{T}, psi::TTN, opstr::String; + kwargs...) where T <: Union{ComplexF64, Float64} N = numsites(psi) sites = get(kwargs, :sites, 1:N) - return map(pos -> expectR(psi, opstr, pos), sites) + return map(pos -> measure(T, psi, opstr, pos), sites) end +measure(psi::TTN, opstr::String; kwargs...) = + measure(ComplexF64, psi, opstr; kwargs) + ################################################################################# -function expectC(psi::TTN, optens::Vector{ITensor})::ComplexF64 +""" + function measure(::Type{ElT} = ComplexF64, psi::TTN, optens::Vector{ITensor}) + +Returns (complex / real) multi-site expectation value (`::ComplexF64` / `::Float64`) +of a given TTN `psi::TTN` for a given vector of single-site operators +(`optens::Vector{ITensor}`). + +For `ElT = Float64`, if the expectation value is complex, raises a warning! +""" +function measure(::Type{T}, psi::TTN, + optens::Vector{ITensor})::T where T <: Union{ComplexF64, Float64} sitenodes = Set{Int2}() for o in optens @@ -67,25 +104,52 @@ function expectC(psi::TTN, optens::Vector{ITensor})::ComplexF64 oc = find_eccentric_central_node(psi.graph, sitenodes) isometrize!(psi, oc) - lnt = LinkTensorsTTN(psi, optens) - phi = psi[oc] - return complex(scalar(product(lnt, psi, phi) * dag(phi))) -end - -################################################################################# + removeqn = false + for o in optens + !hasqns(o) && hasqns(psi) && (removeqn = true) + end + conditional_removeqns(A) = removeqn && hasqns(A) ? removeqns(A) : A -function expectR(psi::TTN, optens::Vector{ITensor})::Float64 - val = expectC(psi, optens) - if abs(imag(val)) > 100 * Float64_threshold() - @warn "`expectR()`: Got complex number with |imaginary part| > $(100 * Float64_threshold()), use `expectC instead !!" + lnt = LinkTensorsTTN(psi, optens) + phi = conditional_removeqns(psi[oc]) + val = complex(scalar(product(lnt, psi, phi) * dag(phi))) + + if T <: Complex + return val + else + if abs(imag(val)) > 100 * Float64_threshold() + @warn "`measure(::Float64)`: Got complex number with |imaginary part| > $(100 * Float64_threshold()) !!" + end + return real(val) end - return real(val) + + return complex(val) end +measure(psi::TTN, optens::Vector{ITensor}) = measure(ComplexF64, psi, optens) + ################################################################################# -function expectC(psi::TTN, oppairs::Vector{Pair{String, Int}}; - isfermions::Bool = true)::ComplexF64 +""" + function measure(::Type{ElT} = ComplexF64, psi::TTN, oppairs::Vector{Pair{String, Int}}; + isfermions::Bool = true) + +Returns (complex / real) multi-site expectation value (`::ComplexF64` / `::Float64`) of +a given TTN `psi::TTN`. `oppairs::Vector{Pair{String, Int}}` contains pairs of operator +names (`String`) and site locations (`Int`). +E.g., for <ψ|OᵢOⱼOₖ... |ψ>, `oppairs = ["O" => i, "O" => j, "O" => k,...]`. + +Fermionic JW strings are added automatically for fermionic operators if +`isfermions::Bool = true` (default). + +For `ElT = Float64`, if the expectation value is complex, raises a warning! + +#### Example: + valueC = measure(psi, ["Cdag" => 2, "C" => 6, "Cdag" => 9, "C" => 12]) + valueR = measure(Float64, psi, ["Cdag" => 2, "C" => 6, "Cdag" => 9, "C" => 12]) +""" +function measure(::Type{T}, psi::TTN, oppairs::Vector{Pair{String, Int}}; + isfermions::Bool = true)::T where T <: Union{ComplexF64, Float64} sites = siteinds(psi) if isfermions @@ -93,18 +157,10 @@ function expectC(psi::TTN, oppairs::Vector{Pair{String, Int}}; end optens = ITensor[op(x.first, siteind(psi, x.second)) for x in oppairs] - return isfermions ? perm * expectC(psi, optens) : expectC(psi, optens) + return isfermions ? perm * measure(T, psi, optens) : measure(T, psi, optens) end +measure(psi::TTN, oppairs::Vector{Pair{String, Int}}; isfermions::Bool = true) = + measure(ComplexF64, psi, oppairs; isfermions) ################################################################################# -function expectR(psi::TTN, oppairs::Vector{Pair{String, Int}}; - isfermions::Bool = true)::Float64 - val = expectC(psi, oppairs; isfermions) - if abs(imag(val)) > 100 * _Float64_Threshold - @warn "`expectR()`: Got complex number with |imaginary part| > $(100 * Float64_threshold()), use `expectC instead !!" - end - return real(val) -end - -################################################################################# diff --git a/src/ttn/optimize_ttn.jl b/src/ttn/optimize_ttn.jl index cf891cc..b83c54c 100644 --- a/src/ttn/optimize_ttn.jl +++ b/src/ttn/optimize_ttn.jl @@ -50,7 +50,7 @@ Base.copy(params::OptimizeParamsTTN) = OptimizeParamsTTN(Base.copy(params.maxdim """ function OptimizeParamsTTN(;maxdim::Vector{Int}, nsweeps::Vector{Int}, - cutoff::Union{Vector{Float64}, Float64} = _Float64_Threshold, + cutoff::Union{Vector{Float64}, Float64} = Float64_threshold(), noise::Union{Vector{Float64}, Float64, Int} = 0.0, noisedecay::Union{Vector{Float64}, Float64, Int} = 1.0, disable_noise_after::Union{Vector{Int}, Int} = typemax(Int)) @@ -71,7 +71,7 @@ Constructor for `OptimizeParamsTTN`. Takes named arguments. throughout the optimization. """ function OptimizeParamsTTN(;maxdim::Vector{Int}, nsweeps::Vector{Int}, - cutoff::Union{Vector{Float64}, Float64} = _Float64_Threshold, + cutoff::Union{Vector{Float64}, Float64} = Float64_threshold(), noise::Union{Vector{Float64}, Float64, Int} = 0.0, noisedecay::Union{Vector{Float64}, Float64, Int} = 1.0, disable_noise_after::Union{Vector{Int}, Int} = typemax(Int)) @@ -109,7 +109,7 @@ end function optimize!(sysenv::StateEnvsTTN, params::OptimizeParamsTTN, sweeppath::Vector{Int2}; - kwargs...)::SweepDataTTN + kwargs...) Performs optimization of the TTN. @@ -143,7 +143,7 @@ See documentation of KrylovKit.jl. - `solver_check_convergence::Bool = false`. #### Return values: - - `SweepData` + - `SweepDataTTN` """ function optimize!(sysenv::StateEnvsTTN, params::OptimizeParamsTTN, @@ -228,17 +228,18 @@ end function optimize(psi0::TTN, H::CouplingModel, params::OptimizeParamsTTN, sweeppath::Vector{Int2}; - kwargs...)::Tuple{Float64, TTN} + kwargs...) function optimize(psi0::TTN, H::CouplingModel, Ms::Vector{TTN}, params::OptimizeParamsTTN, sweeppath::Vector{Int2}; - kwargs...)::Tuple{Float64, TTN} + kwargs...) Performs optimization of the TTN. #### Arguments: - - `sysenv::StateEnvsTTN`. + - `psi0::TTN`: Initial TTN. + - `H::CouplingModel`, `Ms::Vector{MPS}`. - `params::OptimizationParamsTTN`. - `sweeppath::Vector{Int2}`: The path to be followed during optimization sweep. A vector that must contain all the nodes atleast once. diff --git a/src/ttn/sweep_ttn.jl b/src/ttn/sweep_ttn.jl index 60a9f7a..63ef8c4 100644 --- a/src/ttn/sweep_ttn.jl +++ b/src/ttn/sweep_ttn.jl @@ -1,6 +1,21 @@ ################################################################################# +""" + mutable struct SweepDataTTN + sweepcount::Int + maxchi::Vector{Int} + energy::Vector{Float64} + end + +Holds historical data after each (full)sweep of the TTN. Requires for convergence check etc. + - `sweepcount::Int`: Number of fullsweeps performed. + - `maxchi::Vector{Int}`: Maximum MPS bond/link dimensions after every sweep. + - `energy::Vector{Float64}`: Energies after every sweep. + +#### Default constructor: + - `SweepDataTTN()`: Initialize an empty `SweepData` object. +""" mutable struct SweepDataTTN sweepcount::Int maxchi::Vector{Int} @@ -11,6 +26,26 @@ SweepDataTTN() = SweepDataTTN(0, Int[], Float64[]) ################################################################################# +""" + Base.copy(swdata::SweepDataTTN) + +Shallow copy of `SweepDataTTN`. +""" +Base.copy(swdata::SweepDataTTN) = SweepDataTTN(swdata.sweepcount, + Base.copy(swdata.energy), + Base.copy(swdata.entropy) + ) +################################################################################# + +""" + function default_sweeppath(psi::TTN) + +For a given a TTN having the default hierarchical binary tree structure, teturns the default +path for sweeping through the TTN. + +#### Return values: + - `Vector{Int2}`: The list of nodes as the path to be followed during a (half)sweep. +""" function default_sweeppath(psi::TTN) sweeppath = Vector{Int2}() @@ -31,6 +66,54 @@ end ################################################################################# +""" + function fullsweep!(sysenv::StateEnvsTTN, sweeppath::Vector{Int2}, solver, + swdata::SweepDataTTN; kwargs...) + +Perform a fullsweep of the TTN by `solver`. + +#### Arguments: + - `sysenv::StateEnvsTTN`. + - `sweeppath::Vector{Int2}`: The list of nodes as the path to be followed during a (half)sweep. + Must have each node atleast once. For the next halfsweep the reverse path is followed. + - `solver`: Solver for update. Currently only `eig_solver` is supported. + - `swdata::SweepDataTTN`. + +#### Named arguments and their default values: + - `time_step::Union{Float64, ComplexF64, Nothing} = nothing`: Time step for future + functionality. + - `normalize::Bool = true`: Whether to normalize after update. + - `maxdim::Int = typemax(Int)`: Maximum bond dimension after SVD truncation. + - `mindim::Int = 1`: Minimum bond dimension after SVD truncation. + - `cutoff::Float64 = Float64_threshold()`: Cutoff for SVD truncation. + - `svd_alg::String = "divide_and_conquer"`. + - `noise::Float64 = 0.0`. + - `expand_dim::Int = 0` if `noise == 0` else `20`. Dimension to be expanded (on top of `maxdim`) + during subspace expansion. + - `max_expand_dim::Int = 2 * expand_dim`. maximum dimension to be expanded during subspace + expansion. + - `expand_numiter::Int = 4`. Number of iteration for subspace expansion. Should be greater + than 1. + - `linkwise_maxdim::Union{Nothing, Dict{LinkTypeTTN, Int}} = nothing`. Specifies maximum bond + dimension of a particular link. + - `outputlevel::Int = 1`. If `0` prints no information, for `1` outputs after + every fullsweep, if `2` prints at every update step. + +#### Named arguments for `solver` and their default values: +See the documentation of KrylovKit.jl. + - `ishermitian::Bool = true`. + - `solver_tol::Float64 = 1E-14`. + - `solver_krylovdim::Int = 5`. + - `solver_maxiter::Int = 2`. + - `solver_outputlevel::Int = 0`: See `verbosity` in KrylovKit.jl. + - `solver_eager::Bool = false`. + - `solver_check_convergence::Bool = false`. + +#### Return values: + - `::Float64`: Change in Energy ΔE + +`swdata::SweepDataTTN` gets updated. +""" function fullsweep!(sysenv::StateEnvsTTN, sweeppath::Vector{Int2}, solver, swdata::SweepDataTTN; kwargs...) diff --git a/src/ttn/typedef.jl b/src/ttn/typedef.jl index c91b227..c8554cc 100644 --- a/src/ttn/typedef.jl +++ b/src/ttn/typedef.jl @@ -3,7 +3,7 @@ """ const Int2 = Tuple{Int, Int} -Tuple of two `Int`. +Tuple of two `Int`s. """ const Int2 = Tuple{Int, Int} Int2() = (typemin(Int), typemin(Int)) @@ -16,7 +16,7 @@ Int2() = (typemin(Int), typemin(Int)) second::Int2 end -Link type for TTN. +Link / bond between two nodes `first` and `second`. """ struct LinkTypeTTN first::Int2 diff --git a/src/ttn/update_site_ttn.jl b/src/ttn/update_site_ttn.jl index 7c1ac61..cf39029 100644 --- a/src/ttn/update_site_ttn.jl +++ b/src/ttn/update_site_ttn.jl @@ -1,6 +1,45 @@ ################################################################################# +""" + function update_position!(sysenv::StateEnvsTTN, solver, node::Int2; + time_step::Union{Float64, ComplexF64, Nothing}, + normalize::Bool, + maxdim::Int, + mindim::Int, + cutoff::Float64, + svd_alg::String, + kwargs...) + +Moves the orthogonality center to `pos` and update StateEnvsTTN at position `pos` by `solver`. + +#### Arguments: + - `sysenv::StateEnvsTTN` + - `solver`: Solver for update. Currently only `eig_solver` is supported. + - `pos::Int`: Position of the node to be updated. + +#### Named arguments and their default values: + - `time_step::Union{Float64, ComplexF64, Nothing} = nothing`: Time step for future + functionality. + - `normalize::Bool = true`: Whether to normalize after update. + - `maxdim::Int = typemax(Int)`: Maximum bond dimension after SVD truncation. + - `mindim::Int = 1`: Minimum bond dimension after SVD truncation. + - `cutoff::Float64 = Float64_threshold()`: Cutoff for SVD truncation. + - `svd_alg::String = "divide_and_conquer"`. + +#### Named arguments for `solver` and their default values: +See the documentation of KrylovKit.jl. + - `ishermitian::Bool = true`. + - `solver_tol::Float64 = 1E-14`. + - `solver_krylovdim::Int = 5`. + - `solver_maxiter::Int = 2`. + - `solver_outputlevel::Int = 0`: See `verbosity` in KrylovKit.jl. + - `solver_eager::Bool = false`. + - `solver_check_convergence::Bool = false`. + +#### Return values: + - `::Float64`: Energy. +""" function update_position!(sysenv::StateEnvsTTN, solver, node::Int2; time_step::Union{Float64, ComplexF64, Nothing}, normalize::Bool, @@ -25,6 +64,14 @@ end ################################################################################# +""" + function subspace_expand!(psi::TTN, node::Int2, nextnode::Int2, + max_expand_dim::Int, noise::Float64) + +Enlarges the bond domension between `node` and `nextnode` by `max_expand_dim` using the +subspace expansion. +The parameter `noise` controls stength of the perturbation. +""" function subspace_expand!(psi::TTN, node::Int2, nextnode::Int2, max_expand_dim::Int, noise::Float64)::Nothing diff --git a/test/test_MPS_DMRG.jl b/test/test_MPS_DMRG.jl index f677579..9ea793c 100644 --- a/test/test_MPS_DMRG.jl +++ b/test/test_MPS_DMRG.jl @@ -132,6 +132,10 @@ let @show measure(psi_gr, "Sz") @show measure(psi_gr, ["Sz" => 1, "Sz" => 10]) + s1 = removeqns(sites[1]) + s10 = removeqns(sites[10]) + @show measure(psi_gr, [op("Sz", s1), op("Sz", s10)]) + @show measure(psi_gr, [op("Sx", s1), op("Sx", s10)]) @time en1, psi1 = dmrg_ex_2(sites, H, psi0, psi_gr; nsite = nsite) diff --git a/test/test_TTN.jl b/test/test_TTN.jl index c7a5d30..22b63b5 100644 --- a/test/test_TTN.jl +++ b/test/test_TTN.jl @@ -25,7 +25,7 @@ function do_ttn_optimize(sites, H, psi0) sweeppath = default_sweeppath(psi0) - params = OptimizeParamsTTN(; maxdim = [24], nsweeps = [5], + params = OptimizeParamsTTN(; maxdim = [128], nsweeps = [5], cutoff = 1e-14, noise = 1e-2, noisedecay = 5, disable_noise_after = 2) @@ -52,7 +52,13 @@ let sites, H, psi0 = coupling_model() en, psi_gr = do_ttn_optimize(sites, H, psi0) - + @show measure(psi_gr, "Sz") + @show measure(psi_gr, ["Sz" => 1, "Sz" => 10]) + s1 = removeqns(sites[1]) + s10 = removeqns(sites[10]) + @show measure(psi_gr, [op("Sz", s1), op("Sz", s10)]) + @show measure(psi_gr, [op("Sx", s1), op("Sx", s10)]) + en1, psi1 = do_ttn_optimize_ex(sites, H, psi0, psi_gr) end