diff --git a/Project.toml b/Project.toml index 39ab586afe..789a90cc19 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensors" uuid = "9136182c-28ba-11e9-034c-db9fb085ebd5" authors = ["Matthew Fishman ", "Miles Stoudenmire "] -version = "0.8.2" +version = "0.8.3" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/docs/Project.toml b/docs/Project.toml index 763805b88e..daa49d262d 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,7 +1,6 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" -ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67" diff --git a/docs/settings.jl b/docs/settings.jl index 0514c1b62d..69115cdeb0 100644 --- a/docs/settings.jl +++ b/docs/settings.jl @@ -1,18 +1,13 @@ using Documenter using ITensors -using ITensorMPS # Allows using ITensorMPS.jl docstrings in ITensors.jl documentation: # https://github.com/JuliaDocs/Documenter.jl/issues/1734 DocMeta.setdocmeta!(ITensors, :DocTestSetup, :(using ITensors); recursive=true) -DocMeta.setdocmeta!(ITensorMPS, :DocTestSetup, :(using ITensorMPS); recursive=true) sitename = "ITensors.jl" settings = Dict( - # Allows using ITensorMPS.jl docstrings in ITensors.jl documentation: - # https://github.com/JuliaDocs/Documenter.jl/issues/1734 - :modules => [ITensors, ITensorMPS], :pages => [ "Introduction" => "index.md", "Getting Started with ITensor" => [ @@ -21,56 +16,22 @@ settings = Dict( "Enabling Debug Checks" => "getting_started/DebugChecks.md", "Next Steps" => "getting_started/NextSteps.md", ], - "Tutorials" => [ - "DMRG" => "tutorials/DMRG.md", - "Quantum Number Conserving DMRG" => "tutorials/QN_DMRG.md", - "MPS Time Evolution" => "tutorials/MPSTimeEvolution.md", - ], - "Code Examples" => [ - "ITensor Examples" => "examples/ITensor.md", - "MPS and MPO Examples" => "examples/MPSandMPO.md", - "DMRG Examples" => "examples/DMRG.md", - "Physics (SiteType) System Examples" => "examples/Physics.md", - ], - "Documentation" => [ - "Index" => "IndexType.md", - "Index collections" => "IndexSetType.md", - "ITensor" => "ITensorType.md", - "MPS and MPO" => "MPSandMPO.md", - "QN" => "QN.md", - "SiteType and op, state, val functions" => "SiteType.md", - "SiteTypes Included with ITensor" => "IncludedSiteTypes.md", - "DMRG" => [ - "DMRG.md", - "Sweeps.md", - "ProjMPO.md", - "ProjMPOSum.md", - "Observer.md", - "DMRGObserver.md", - ], - "OpSum" => "OpSum.md", - ], + "Code Examples" => ["ITensor Examples" => "examples/ITensor.md"], + "Documentation" => + ["Index" => "IndexType.md", "ITensor" => "ITensorType.md", "QN" => "QN.md"], "Frequently Asked Questions" => [ "Programming Language (Julia, C++, ...) FAQs" => "faq/JuliaAndCpp.md", - "DMRG FAQs" => "faq/DMRG.md", - "Quantum Number (QN) FAQs" => "faq/QN.md", "ITensor Development FAQs" => "faq/Development.md", - "Relationship of ITensor to other tensor libraries FAQs" => "faq/RelationshipToOtherLibraries.md", "Julia Package Manager FAQs" => "faq/JuliaPkg.md", "High-Performance Computing FAQs" => "faq/HPC.md", ], "Upgrade guides" => ["Upgrading from 0.1 to 0.2" => "UpgradeGuide_0.1_to_0.2.md"], - "ITensor indices and Einstein notation" => "Einsum.md", "Advanced Usage Guide" => [ - "Advanced Usage Guide" => "AdvancedUsageGuide.md", "Multithreading" => "Multithreading.md", "Running on GPUs" => "RunningOnGPUs.md", - "Symmetric (QN conserving) tensors: background and usage" => "QNTricks.md", - "Timing and profiling" => "CodeTiming.md", "Contraction sequence optimization" => "ContractionSequenceOptimization.md", "HDF5 File Formats" => "HDF5FileFormats.md", ], - "Developer Guide" => "DeveloperGuide.md", ], :format => Documenter.HTML(; assets=["assets/favicon.ico"], prettyurls=false), :doctest => true, diff --git a/docs/src/AdvancedUsageGuide.md b/docs/src/AdvancedUsageGuide.md deleted file mode 100644 index 54d3da6660..0000000000 --- a/docs/src/AdvancedUsageGuide.md +++ /dev/null @@ -1,1124 +0,0 @@ -# [Advanced ITensor Usage Guide](@id advanced_usage_guide) - -## Installing and updating ITensors.jl - -The ITensors package can be installed with the Julia package manager. -Assuming you have already downloaded Julia, which you can get -[here](https://julialang.org/downloads/), from the Julia REPL, -type `]` to enter the Pkg REPL mode and run: -``` -$ julia -``` -```julia -julia> ] - -pkg> add ITensors -``` - -Or, equivalently, via the `Pkg` API: - -```julia -julia> import Pkg; Pkg.add("ITensors") -``` - -We recommend using ITensors.jl with Intel MKL in order to get the -best possible performance. If you have not done so already, you can -replace the current BLAS and LAPACK implementation used by Julia with -MKL by using the MKL.jl package. Please follow the instructions -[here](https://github.com/JuliaComputing/MKL.jl). - -To use the latest registered (stable) version of ITensors.jl, use `update ITensors` -in `Pkg` mode or `import Pkg; Pkg.update("ITensors")`. -We will commonly release new patch versions (such as updating from `v0.1.12` to -`v0.1.13`) with bug fixes and improvements. However, make sure to double check before -updating between minor versions (such as from `v0.1.41` to `v0.2.0`) because new minor -releases may be breaking. - -Remember that if you are compiling system images of ITensors.jl, such as with the -`ITensors.compile()` command, you will need to rerurn this command to compile -the new version of ITensor after an update. - -To try the "development branch" of ITensors.jl (for example, if -there is a feature or fix we added that hasn't been released yet), -you can do `add ITensors#main`. You can switch back to the latest -released version with `add ITensors`. Using the development/main -branch is generally not encouraged unless you know what you are doing. - -## Using ITensors.jl in the REPL - -There are many ways you can write code based on ITensors.jl, ranging -from using it in the REPL to writing a small script to making a -package that depends on it. - -For example, you can just start the REPL from your command line like: -``` -$ julia -``` -assuming you have an available version of Julia with the ITensors.jl -package installed. Then just type: -```julia -julia> using ITensors -``` -and start typing ITensor commands. For example: -```julia -julia> i = Index(2, "i") -(dim=2|id=355|"i") - -julia> A = random_itensor(i, i') -ITensor ord=2 (dim=2|id=355|"i") (dim=2|id=355|"i")' -NDTensors.Dense{Float64,Array{Float64,1}} - -julia> @show A; -A = ITensor ord=2 -Dim 1: (dim=2|id=355|"i") -Dim 2: (dim=2|id=355|"i")' -NDTensors.Dense{Float64,Array{Float64,1}} - 2×2 - 1.2320011464276275 1.8504245734277216 - 1.0763652402177477 0.030353720156277037 - -julia> (A*dag(A))[] -3.9627443142240617 -``` - -Note that there are some "gotchas" with working in the REPL like this. -Technically, all commands in the REPL are in the "global scope". -The global scope might not work as you would expect, for example: -```julia -julia> for _ in 1:3 - A *= 2 - end -ERROR: UndefVarError: A not defined -Stacktrace: - [1] top-level scope at ./REPL[12]:2 - [2] eval(::Module, ::Any) at ./boot.jl:331 - [3] eval_user_input(::Any, ::REPL.REPLBackend) at /home/mfishman/software/julia-1.4.0/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86 - [4] run_backend(::REPL.REPLBackend) at /home/mfishman/.julia/packages/Revise/AMRie/src/Revise.jl:1023 - [5] top-level scope at none:0 -``` -since the `A` inside the for-loop introduces a new local variable. -Some alternatives are to wrap that part of the code in a let-block -or a function: -```julia -julia> function f(A) - for _ in 1:3 - A *= 2 - end - A - end -f (generic function with 1 method) - -julia> A = f(A) -ITensor ord=2 (dim=2|id=355|"i") (dim=2|id=355|"i")' -NDTensors.Dense{Float64,Array{Float64,1}} - -julia> @show A; -A = ITensor ord=2 -Dim 1: (dim=2|id=355|"i") -Dim 2: (dim=2|id=355|"i")' -NDTensors.Dense{Float64,Array{Float64,1}} - 2×2 - 9.85600917142102 14.803396587421773 - 8.610921921741982 0.2428297612502163 -``` -In this particular case, you can alternatively modify the ITensor -in-place: -```julia -julia> for _ in 1:3 - A ./= 2 - end - -julia> @show A; -A = ITensor ord=2 -Dim 1: (dim=2|id=355|"i") -Dim 2: (dim=2|id=355|"i")' -NDTensors.Dense{Float64,Array{Float64,1}} - 2×2 - 1.2320011464276275 1.8504245734277216 - 1.0763652402177477 0.030353720156277037 -``` - -A common place you might accidentally come across this is when -you are creating a Hamiltonian with `OpSum`: -```julia -julia> using ITensors, ITensorMPS - -julia> N = 4; - -julia> sites = siteinds("S=1/2",N); - -julia> os = OpSum(); - -julia> for j=1:N-1 - os += "Sz", j, "Sz", j+1 - end -ERROR: UndefVarError: os not defined -Stacktrace: - [1] top-level scope at ./REPL[16]:2 - [2] eval(::Module, ::Any) at ./boot.jl:331 - [3] eval_user_input(::Any, ::REPL.REPLBackend) at /home/mfishman/software/julia-1.4.0/share/julia/stdlib/v1.4/REPL/src/REPL.jl:86 - [4] run_backend(::REPL.REPLBackend) at /home/mfishman/.julia/packages/Revise/AMRie/src/Revise.jl:1023 - [5] top-level scope at none:0 -``` -In this case, you can use `os .+= ("Sz", j, "Sz", j+1)`, -`add!(os, "Sz", j, "Sz", j+1)`, or wrap your code in a let-block -or function. - -Take a look at Julia's documentation [here](https://docs.julialang.org/en/v1/manual/variables-and-scoping/) -for rules on scoping. Also note that this behavior is particular -to Julia v1.4 and below, and is expected to change in v1.5. - -Note that the REPL is very useful for prototyping code quickly, -but working directly in the REPL and outside of functions can -cause sub-optimal performance. See Julia's [performance tips](https://docs.julialang.org/en/v1/manual/performance-tips/index.html) -for more information. - -We recommend the package [OhMyREPL](https://kristofferc.github.io/OhMyREPL.jl/latest/) which adds syntax highlighting to the Julia REPL. - -## Finding documentation interactively - -Julia provides many tools for searching for documentation interactively at the REPL. Say that you want to learn more about how to use an ITensor from the command line. You can start by typing `?` followed by `ITensor`: -```julia -julia> using ITensors - -julia> ?ITensor -search: ITensor ITensors itensor random_itensor - - An ITensor is a tensor whose interface is independent of its - memory layout. Therefore it is not necessary to know the ordering - of an ITensor's indices, only which indices an ITensor has. - Operations like contraction and addition of ITensors automatically - handle any memory permutations. - - Examples - ≡≡≡≡≡≡≡≡≡≡ - - julia> i = Index(2, "i") - (dim=2|id=287|"i") - - julia> A = random_itensor(i', i) - ITensor ord=2 (dim=2|id=287|"i")' (dim=2|id=287|"i") - NDTensors.Dense{Float64,Array{Float64,1}} - - julia> @show A; - A = ITensor ord=2 - Dim 1: (dim=2|id=287|"i")' - Dim 2: (dim=2|id=287|"i") - NDTensors.Dense{Float64,Array{Float64,1}} - 2×2 - 0.28358594718392427 1.4342219756446355 - 1.6620103556283987 -0.40952231269251566 - - julia> @show inds(A); - inds(A) = IndexSet{2} (dim=2|id=287|"i")' (dim=2|id=287|"i") -[...] -``` -(the specific output may be different for different versions of ITensors.jl as we update the docs). You can use the help prompt (which you get by typing `?` at the REPL) to print out documentation for types and methods. - -Another way to get information about types is with the function `fieldnames`: -```julia -julia> fieldnames(ITensor) -(:store, :inds) -``` -which shows the fields of a type. Note that in general the specific names of the fields and structures of types may change (we consider those to be internal details), however we often make functions to access the fields of a type that have the same name as the field, so it is a good place to get started. For example, you can access the storage and indices of an ITensor `A` with the functions `store(A)` and `inds(A)`. - -Another helpful function is `apropos`, which search through all documentation for a string (ignoring the case) and prints a list of all types and methods with documentation that contain the string. - -Based on the `apropos` function, we can make some helper functions that may be useful. For example: -```julia -using ITensors - -function finddocs(s) - io = IOBuffer() - apropos(io, s) - v = chomp(String(take!(io))) - return split(v, "\n") -end - -function finddocs(s...) - intersect(finddocs.(s)...) -end - -found_methods = finddocs("indices", "set difference") -display(found_methods) -``` -returns: -```julia -3-element Array{SubString{String},1}: - "ITensors.noncommoninds" - "Base.setdiff" - "ITensors.uniqueinds" -``` -which are the functions that have docs that contain the strings `"indices"` and `"set difference"`. We can print the docs for `uniqueinds` to find: -```julia -help?> uniqueinds -search: uniqueinds unique_siteinds uniqueind uniqueindex - - uniqueinds(A, B; kwargs...) - uniqueinds(::Order{N}, A, B; kwargs...) - - - Return an IndexSet with indices that are unique to the set of - indices of A and not in B (the set difference). - - Optionally, specify the desired number of indices as Order(N), - which adds a check and can be a bit more efficient. -``` - -We can also filter the results to only specify functions from certain modules, for example: -```julia -julia> filter(x -> startswith(x, "ITensors"), finddocs("indices", "set difference")) -2-element Array{SubString{String},1}: - "ITensors.noncommoninds" - "ITensors.uniqueinds" - -julia> filter(x -> !startswith(x, "ITensors"), finddocs("indices", "set difference")) -1-element Array{SubString{String},1}: - "Base.setdiff" -``` -Ideally we could have `apropos` do a "smart" Google-like search of the appropriate docstrings, but this is a pretty good start. - -Additionally, the `names` function can be useful, which prints the names of all functions and types that are exported by a module. For example: -```julia -julia> names(ITensors) -264-element Array{Symbol,1}: - Symbol("@OpName_str") - Symbol("@SiteType_str") - Symbol("@StateName_str") - Symbol("@TagType_str") - Symbol("@disable_warn_order") - Symbol("@reset_warn_order") - Symbol("@set_warn_order") - Symbol("@ts_str") - :AbstractObserver - :OpSum - :DMRGObserver - :ITensor - :ITensors - :Index -[...] -``` -Of course this is a very long list (and the methods are returned as `Symbol`s, which are like strings but not as easy to work with). However, we can convert the list to strings and filter the strings to find functions we are interested in, for example: -```julia -julia> filter(x -> contains(x, "common") && contains(x, "ind"), String.(names(ITensors))) -8-element Array{String,1}: - "common_siteind" - "common_siteinds" - "commonind" - "commonindex" - "commoninds" - "hascommoninds" - "noncommonind" - "noncommoninds" -``` - -Julia types do not have member functions, so people coming from object oriented programming languages may find that at first it is more difficult to find methods that are applicable to a certain type. However, Julia has many fantastic tools for introspection that we can use to make this task easier. - -## Make a small project based on ITensors.jl - -Once you start to have longer code, you will want to put your -code into one or more files. For example, you may have a short script -with one or more functions based on ITensors.jl: -```julia -# my_itensor_script.jl -using ITensors - -function norm2(A::ITensor) - return (A*dag(A))[] -end -``` -Then, in the same directory as your script `my_itensor_script.jl`, -just type: -```julia -julia> include("my_itensor_script.jl"); - -julia> i = Index(2; tags="i"); - -julia> A = random_itensor(i', i); - -julia> norm2(A) -[...] -``` - -As your code gets longer, you can split it into multiple files -and `include` this files into one main project file, for example -if you have two files with functions in them: -```julia -# file1.jl - -function norm2(A::ITensor) - return (A*dag(A))[] -end -``` -and -```julia -# file2.jl - -function square(A::ITensor) - return A .^ 2 -end -``` - -```julia -# my_itensor_project.jl - -using ITensors - -include("file1.jl") - -include("file2.jl") -``` -Then, as before, you can use your functions at the Julia REPL -by just including the file `my_itensor_project.jl`: -```julia -julia> include("my_itensor_project.jl"); - -julia> i = Index(2; tags="i"); - -julia> A = random_itensor(i', i); - -julia> norm2(A) -[...] - -julia> square(A) -[...] -``` - -As your code gets more complicated and has more files, it is helpful -to organize it into a package. That will be covered in the -next section. - -## Make a Julia package based on ITensors.jl - -In this section, we will describe how to make a Julia package based on -ITensors.jl. This is useful to do when your project gets longer, -since it helps with: - - Code organization. - - Adding dependencies that will get automatically installed through Julia's package system. - - Versioning. - - Automated testing. - - Code sharing and easier package installation. - - Officially registering your package with Julia. -and many more features that we will mention later. - -Start up Julia and install [PkgTemplates](https://invenia.github.io/PkgTemplates.jl/stable/) -```julia -$ julia - -julia> ] - -pkg> add PkgTemplates -``` -then press backspace and type: -``` -julia> using PkgTemplates - -julia> t = Template(; user="your_github_username", plugins=[Git(; ssh=true),]) - -julia> t("MyITensorsPkg") -``` -You should put your Github account name instead of `"your_github_username"`, -if you want to use Github to host your package. -The option `plugins=[Git(; ssh=true),]` sets the Github authentication to use -ssh, which is generally more convenient. You can switch to https (where you -have to type your username and password to push changes) by setting `ssh=false` -or leaving off `plugins=[...]`. By default, the package will be located in -the directory `~/.julia/dev`, you can change this with the keyword argument -`dir=[...]`. However, `~/.julia/dev` is recommended since that is the directory -Julia's package manager (and other packages like `Revise`) will look for development -packages. Please see the `PkgTemplate` documentation for more customization options. - -Then, we want to tell Julia about our new package. We do this as -follows: -```julia -julia> ] - -pkg> dev ~/.julia/dev/MyITensorsPkg -``` -then you can do: -```julia -julia> using MyITensorsPkg -``` -from any directory to use your new package. However, it doesn't -have any functions available yet. Additionally, there should be -an empty test file already set up here: -``` -~/.julia/dev/MyITensorsPkg/test/runtests.jl -``` -which you can run from any directory like: -```julia -julia> ] - -pkg> test MyITensorsPkg -``` -It should show something like: -```julia -[...] -Test Summary: | -MyITensorsPkg.jl | No tests - Testing MyITensorsPkg tests passed -``` -since there are no tests yet. - -First we want to add ITensors as a dependency of our package. -We do this by "activating" our package environment and then -adding ITensors: -```julia -julia> ] - -pkg> activate MyITensorsPkg - -(MyITensorsPkg) pkg> add ITensors -``` -This will edit the file `~/.julia/dev/MyITensorsPkg/Project.toml` -and add the line -``` -[deps] -ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" -``` -Because your package is under development, back in the main -Pkg environment you should type `resolve`: -```julia -(MyITensorsPkg) pkg> activate - -pkg> resolve -``` -Now, if you or someone else uses the package, it will automatically -install ITensors.jl for you. - -Now your package is set up to develop! Try editing the file -`~/.julia/dev/MyITensorsPkg/src/MyITensorsPkg.jl` and add the -`norm2` function, which calculates the squared norm of an ITensor: -```julia -module MyITensorsPkg - -using ITensors - -export norm2 - -norm2(A::ITensor) = (A*dag(A))[] - -end -``` -The export command makes `norm2` available in the namespace without -needing to type `MyITensorsPkg.norm2` when you do -`using MyITensorsPkg`. Now in a new Julia session you can do: -```julia -julia> using ITensors - -julia> i = Index(2) -(dim=2|id=263) - -julia> A = random_itensor(i) -ITensor ord=1 (dim=2|id=263) -NDTensors.Dense{Float64,Array{Float64,1}} - -julia> norm(A)^2 -6.884457016011188 - -julia> norm2(A) -ERROR: UndefVarError: norm2 not defined -[...] - -julia> using MyITensorsPkg - -julia> norm2(A) -6.884457016011188 -``` -Unfortunately, if you continue to edit the file `MyITensorsPkg.jl`, -even if you type `using MyITensorsPkg` again, if you are in the -same Julia session the changes will not be reflected, and -you will have to restart your Julia session. The -[Revise](https://timholy.github.io/Revise.jl/stable/) package -will allow you to edit your package files and have the changes -reflected in real time in your current Julia session, so you -don't have to restart the session. - -Now, we can add some tests for our new functionality. Edit -the file `~/.julia/dev/MyITensorsPkg/test/runtests.jl` to -look like: -```julia -using MyITensorsPkg -using ITensors -using Test - -@testset "MyITensorsPkg.jl" begin - i = Index(2) - A = random_itensor(i) - @test isapprox(norm2(A), norm(A)^2) -end -``` -Now when you test your package you should see: -```julia -pkg> test MyITensorsPkg -[...] -Test Summary: | Pass Total -MyITensorsPkg.jl | 1 1 - Testing MyITensorsPkg tests passed -``` - -Your package should already be set up as a git repository by -the `PkgTemplates` commands we started with. -We recommend using Github or similar versions control systems -for your packages, especially if you plan to make them public -and officially register them as Julia packages. - -You can set up your local package as a Github repository by -following the steps [here](https://help.github.com/en/github/importing-your-projects-to-github/adding-an-existing-project-to-github-using-the-command-line). Many of the steps may be unnecessary since they -were already set up by `PkgTemplates`. You should be able to -go to the website [here](https://github.com/new), create a new -Github repository with the name `MyITensorsPkg.jl`, and then following -the instructions under "push an existing repository from the command line". - -You may also want to switch between HTTPS and SSH authentication -as described [here](https://help.github.com/en/github/using-git/changing-a-remotes-url), -if you didn't choose your preferred authentication protocol with -PkgTemplates. - -There are many more features you can add to your package through -various Julia packages and Github, for example: - - Control of precompilation with tools like [SnoopCompile](https://timholy.github.io/SnoopCompile.jl/stable/). - - Automatic testing of your package at every pull request/commit with Github Actions, Travis, or similar services. - - Automated benchmarking of your package at every pull request with [BenchmarkTools](https://github.com/JuliaCI/BenchmarkTools.jl), [PkgBenchmark](https://juliaci.github.io/PkgBenchmark.jl/stable/) and [BenchmarkCI](https://github.com/tkf/BenchmarkCI.jl). - - Automated building of your documentation with [Documenter](https://juliadocs.github.io/Documenter.jl/stable/). - - Compiling your package with [PackageCompiler](https://julialang.github.io/PackageCompiler.jl/dev/). - - Automatically check what parts of your code your tests check with code coverage. - - Officially register your Julia package so that others can easily install it and follow along with updated versions using the [Registrator](https://juliaregistries.github.io/Registrator.jl/stable/). -You can take a look at the [ITensors](https://github.com/ITensor/ITensors.jl) -Github page for inspiration on setting up some of these services -and ideas for organizing your package. - -## Developing ITensors.jl - -This section is for someone who is interested in modifying the source code of ITensors.jl, -and then possibly contribute you changes to the official ITensors.jl package. - -This should not be necessary for most people. If for whatever reason you think that the functionality -of ITensors.jl needs to be modified, oftentimes you can add new functions outside of ITensors.jl -or directly overload a function of ITensors.jl (for example with the [import](https://docs.julialang.org/en/v1/manual/modules/#using-and-import-with-specific-identifiers,-and-adding-methods) keyword). - -However, if you would like to only modify parts of the internals of an ITensors.jl function, -and/or plan to contribute changes like bug fixes or new features to the official ITensors.jl -package, this section is for you. - -If you install a package like ITensors with the package manager using the standard `Pkg.add` command: -```julia -julia> using Pkg - -julia> Pkg.add("ITensors") -``` -it will automatically clone the latest registered/tagged version of `ITensors` in a randomly -generated directory inside `~/.julia/packages`. You can find out what version you are using with `Pkg.status`: -```julia -julia> Pkg.status("ITensors") - Status `~/.julia/environments/v1.7/Project.toml` - [9136182c] ITensors v0.2.16 -``` -and you can use [`pkgdir`](https://docs.julialang.org/en/v1/base/base/#Base.pkgdir-Tuple{Module}) -to find out the directory of the source code of a package that you have loaded: -```julia -julia> using ITensors - -julia> pkgdir(ITensors) -"/home/mfishman/.julia/packages/ITensors/cu9Bo" -``` -The source code of a package loaded in this way is read-only, so you won't be able to modify it. - -If you want to modify the source code of `ITensors.jl`, you should check out the packages -`NDTensors.jl` and `ITensors.jl` in development mode with `Pkg.develop`: -```julia -julia> Pkg.develop(["NDTensors", "ITensors"]) -Path `/home/mfishman/.julia/dev/ITensors` exists and looks like the correct repo. Using existing path. - Resolving package versions... - Updating `~/.julia/environments/v1.7/Project.toml` - [9136182c] ~ ITensors v0.2.16 ⇒ v0.2.16 `~/.julia/dev/ITensors` - [23ae76d9] ~ NDTensors v0.1.35 ⇒ v0.1.35 `~/.julia/dev/ITensors/NDTensors` - Updating `~/.julia/environments/v1.7/Manifest.toml` - [9136182c] ~ ITensors v0.2.16 ⇒ v0.2.16 `~/.julia/dev/ITensors` - [23ae76d9] ~ NDTensors v0.1.35 ⇒ v0.1.35 `~/.julia/dev/ITensors/NDTensors` - -julia> Pkg.status(["NDTensors", "ITensors"]) - Status `~/.julia/environments/v1.7/Project.toml` - [9136182c] ITensors v0.2.16 `~/.julia/dev/ITensors` - [23ae76d9] NDTensors v0.1.35 `~/.julia/dev/ITensors/NDTensors` -``` -Then, Julia will use the version of `ITensors.jl` living in the directory `~/.julia/dev/ITensors` -and the version of `NDTensors.jl` living in the directory `~/.julia/dev/ITensors/NDTensors`, -though you may need to restart Julia for this to take affect. - -We recommend checking out the development versions of both `NDTensors.jl` and `ITensors.jl` -since we often develop both packages tandem, so the development branch -of `ITensors.jl` may rely on changes we make in `NDTensors.jl`. - -By default, when you modify code in `~/.julia/dev/ITensors` or `~/.julia/dev/ITensors/NDTensors` -you will need to restart Julia for the changes to take affect. -A way around this issue is the [Revise](https://timholy.github.io/Revise.jl/stable/) package. -We highly recommend using the [Revise](https://timholy.github.io/Revise.jl/stable/) package -when you are developing packages, which automatically detects changes you are making to -a package you have checked out for development and edit code and not have to restart your Julia session. -In short, if you have `Revise.jl` loaded, you can edit the code in `~/.julia/dev/ITensors` -or `~/.julia/dev/ITensors/NDTensors` and the changes you make will be reflected on the fly as -you use the package (there are some limitations, for example you will need to restart Julia -if you change the definitions of types). - -Note that the code in `~/.julia/dev/ITensors` is just a git repository cloned from -the repository https://github.com/ITensor/ITensors.jl, so you can do anything that -you would with any other git repository (use forks of the project, check out branches, -push and pull changes, etc.). - -The standard procedure for submitting a bug fix or new feature to ITensors.jl would then -be to first [fork the ITensors.jl repository](https://docs.github.com/en/get-started/quickstart/fork-a-repo). -Then, check out your fork for development with: -```julia -julia> using Pkg - -julia> Pkg.develop(url="https://github.com/mtfishman/ITensors.jl") -``` -where you would replace `mtfishman` with your own Github username. -Make the changes to the code in `~/.julia/dev/ITensors`, push the changes to your fork, and then -[make a pull request](https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/proposing-changes-to-your-work-with-pull-requests/creating-a-pull-request) to the [ITensors.jl Github repository](https://github.com/ITensor/ITensors.jl/compare). - -To go back to the official version of the `NDTensors.jl` and `ITensors.jl` packages, you can use the command `Pkg.free(["NDTensors", "ITensors"])`: -```julia -julia> Pkg.free(["NDTensors", "ITensors"]) - Resolving package versions... - Updating `~/.julia/environments/v1.7/Project.toml` - [9136182c] ~ ITensors v0.2.16 `~/.julia/dev/ITensors` ⇒ v0.2.16 - [23ae76d9] ~ NDTensors v0.1.35 `~/.julia/dev/ITensors/NDTensors` ⇒ v0.1.35 - Updating `~/.julia/environments/v1.7/Manifest.toml` - [9136182c] ~ ITensors v0.2.16 `~/.julia/dev/ITensors` ⇒ v0.2.16 - [23ae76d9] ~ NDTensors v0.1.35 `~/.julia/dev/ITensors/NDTensors` ⇒ v0.1.35 - -julia> Pkg.status(["NDTensors", "ITensors"]) - Status `~/.julia/environments/v1.7/Project.toml` - [9136182c] ITensors v0.2.16 - [23ae76d9] NDTensors v0.1.35 -``` -so it returns to the version of the package you would have just after installing with `Pkg.add`. - -Some of the Julia package development workflow definitely takes some getting used to, -but once you figure out the "flow" and have a picture of what is going on there are -only a small set of commands you really need to use. - -A small note is that we follow the [Blue style guide](https://github.com/invenia/BlueStyle) -for formatting the source code in ITensors.jl. -To make this more automated, we use the wonderful package -[JuliaFormatter.jl](https://github.com/domluna/JuliaFormatter.jl). -To format your developed version of ITensors.jl, all you have to do is change your directory -to `~/.julia/dev/ITensors` and run the command `format(".")` after loading the `JuliaFormatter` -package: -```julia -julia> using Pkg - -julia> Pkg.status("ITensors") - Status `~/.julia/environments/v1.7/Project.toml` - [9136182c] ITensors v0.2.16 `~/.julia/dev/ITensors` - -julia> using ITensors - -julia> pkgdir(ITensors) -"/home/mfishman/.julia/dev/ITensors" - -julia> cd(pkgdir(ITensors)) - -julia> using JuliaFormatter - -julia> format(".") -false - -julia> format(".") # Check the formatting succeeded -true -``` -This will automatically change the style of the code according to the `Blue` style guide. -The `format` command returns `false` if the code was not already formatted (and therefore -if the command made changes to the source code to follow the style guide), and -returns `true` otherwise. - -If you make changes to ITensors that you think will be useful to others, such as fixing bugs -or adding new features, please consider making a [pull request](https://github.com/ITensor/ITensors.jl/compare). -However, please ask us first before doing so -- either by raising an [issue on Github](https://github.com/ITensor/ITensors.jl/issues) or asking a question on the [ITensor support forum](http://itensor.org/support/) -- -to make sure it is a change or addition that we will want to include or to check that it is not something -we are currently working on. Coordinating with us in that way will help save your time and energy as well as ours! - -[Here](https://www.youtube.com/watch?v=QVmU29rCjaA) is a great introduction to Julia package development -as well as making pull requests to existing Julia packages by the irreplacable Chris Rackauckas. - -## Compiling ITensors.jl - -You might notice that the time to load ITensors.jl (with `using -ITensors`) and the time to run your first few ITensor commands is -slow. This is due to Julia's just-in-time (JIT) compilation. -Julia is compiling special versions of each function that is -being called based on the inputs that it gets at runtime. This -allows it to have fast code, often nearly as fast as fully compiled -languages like C++, while still being a dynamic language. - -However, the long startup time can still be annoying. In this section, -we will discuss some strategies that can be used to minimize this -annoyance, for example: - - Precompilation. - - Staying in the same Julia session with [Revise](https://timholy.github.io/Revise.jl/stable/). - - Using [PackageCompiler](https://julialang.github.io/PackageCompiler.jl/dev/) to compile ITensors.jl ahead of time. - -Precompilation is performed automatically when you first install -ITensors.jl or update a version and run the command `using ITensors` -for the first time. For example, when you first use ITensors after -installation or updating, you will see: -```julia -julia> using ITensors -[ Info: Precompiling ITensors [9136182c-28ba-11e9-034c-db9fb085ebd5] -``` -The process is done automatically, and -puts some compiled binaries in your `~/.julia` directory. The -goal is to decrease the time it takes when you first type -`using ITensors` in your next Julia session, and also the time -it takes for you to first run ITensor functions in a new -Julia session. This helps the startup time, but currently doesn't -help enough. -This is something both ITensors.jl and the Julia language will try -to improve over time. - -To avoid this time, it is recommended that you work as much as you -can in a single Julia session. You should not need to restart your -Julia session very often. For example, if you are writing code in -a script, just `include` the file again which will pull in the new -changes to the script (the exception is if you change the definition -of a type you made, which would requiring restarting the REPL). - -If you are working on a project, we highly recommend using the -[Revise](https://timholy.github.io/Revise.jl/stable/) package -which automatically detects changes you are making in your -packages and reflects them real-time in your current REPL session. -Using these strategies should minimize the number of times you -need to restart your REPL session. - -If you plan to use ITensors.jl directly from the command line -(i.e. not from the REPL), and the startup time is an issue, -you can try compiling ITensors.jl using [PackageCompiler](https://julialang.github.io/PackageCompiler.jl/dev/). - -Before using PackageCompiler to compile ITensors, when we first start using ITensors.jl we might see: -```julia -julia> @time using ITensors - 3.845253 seconds (10.96 M allocations: 618.071 MiB, 3.95% gc time) - -julia> @time i = Index(2); - 0.000684 seconds (23 allocations: 20.328 KiB) - -julia> @time A = random_itensor(i', i); - 0.071022 seconds (183.24 k allocations: 9.715 MiB) - -julia> @time svd(A, i'); - 5.802053 seconds (24.56 M allocations: 1.200 GiB, 7.83% gc time) - -julia> @time svd(A, i'); - 0.000177 seconds (450 allocations: 36.609 KiB) -``` -ITensors provides the command `ITensors.compile()` to create what is -called a "custom system image", a custom version of Julia that -includes a compiled version of ITensors (see the [PackageCompiler documentation](https://julialang.github.io/PackageCompiler.jl/dev/) for more details). - -!!! compat "ITensors 0.7" - As of ITensors 0.7, you must now install and load both the - [ITensorMPS.jl](https://github.com/ITensor/ITensorMPS.jl) package - and the [PackageCompiler.jl](https://github.com/JuliaLang/PackageCompiler.jl) - package in order to use `ITensors.compile()`, since it relies on running MPS/MPO - functionality as example code for Julia to compile and is based in a package - extension in order to make `PackageCompiler.jl` an optional dependency. - -Just run the commands: -``` -julia> using ITensors, ITensorMPS - -julia> using PackageCompiler - -julia> ITensors.compile() -[...] -``` -By default, this will create the file `sys_itensors.so` in the directory -`~/.julia/sysimages`. -Then if we start julia with: -``` -$ julia --sysimage ~/.julia/sysimages/sys_itensors.so -``` -then you should see something like: -```julia -julia> @time using ITensors - 0.330587 seconds (977.61 k allocations: 45.807 MiB, 1.89% gc time) - -julia> @time i = Index(2); - 0.000656 seconds (23 allocations: 20.328 KiB) - -julia> @time A = random_itensor(i', i); - 0.000007 seconds (7 allocations: 576 bytes) - -julia> @time svd(A, i'); - 0.263526 seconds (290.02 k allocations: 14.220 MiB) - -julia> @time svd(A, i'); - 0.000135 seconds (350 allocations: 29.984 KiB) -``` -which is much better. - -Note that you will have to recompile ITensors with the command -`ITensors.compile()` any time that you update the version of ITensors -in order to keep the system image updated. We hope to make this -process more automated in the future. - -## Benchmarking and profiling - -Julia has great built-in tools for benchmarking and profiling. -For benchmarking fast code at the command line, you can use -[BenchmarkTools](https://github.com/JuliaCI/BenchmarkTools.jl/blob/main/doc/manual.md): -```julia -julia> using ITensors; - -julia> using BenchmarkTools; - -julia> i = Index(100, "i"); - -julia> A = random_itensor(i, i'); - -julia> @btime 2*$A; - 4.279 μs (8 allocations: 78.73 KiB) -``` - -We recommend packages like [ProfileView](https://github.com/timholy/ProfileView.jl) -to get detailed profiles of your code, in order to pinpoint functions -or lines of code that are slower than they should be. - -## ITensor type design and writing performant code - -Advanced users might notice something strange about the definition -of the ITensor type, that it is often not "type stable". Some of -this is by design. The definition for ITensor is: -```julia -mutable struct ITensor - inds::IndexSet - store::TensorStorage -end -``` -These are both abstract types, which is something that is generally -discouraged for peformance. - -This has a few disadvantages. Some code that you might expect to be -type stable, like `getindex`, is not, for example: -```julia -julia> i = Index(2, "i"); - -julia> A = random_itensor(i, i'); - -julia> @code_warntype A[i=>1, i'=>2] -Variables - #self#::Core.Compiler.Const(getindex, false) - T::ITensor - ivs::Tuple{Pair{Index{Int64},Int64}} - p::Tuple{Union{Nothing, Int64}} - vals::Tuple{Any} - -Body::Number -1 ─ %1 = NDTensors.getperm::Core.Compiler.Const(NDTensors.getperm, false) -│ %2 = ITensors.inds(T)::IndexSet{1,IndexT,DataT} where DataT<:Tuple where IndexT<:Index -│ %3 = Base.broadcasted(ITensors.ind, ivs)::Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(ind),Tuple{Tuple{Pair{Index{Int64},Int64}}}} -│ %4 = Base.materialize(%3)::Tuple{Index{Int64}} -│ (p = (%1)(%2, %4)) -│ %6 = NDTensors.permute::Core.Compiler.Const(NDTensors.permute, false) -│ %7 = Base.broadcasted(ITensors.val, ivs)::Base.Broadcast.Broadcasted{Base.Broadcast.Style{Tuple},Nothing,typeof(val),Tuple{Tuple{Pair{Index{Int64},Int64}}}} -│ %8 = Base.materialize(%7)::Tuple{Int64} -│ (vals = (%6)(%8, p)) -│ %10 = Core.tuple(T)::Tuple{ITensor} -│ %11 = Core._apply_iterate(Base.iterate, Base.getindex, %10, vals)::Number -│ %12 = Core.typeassert(%11, ITensors.Number)::Number -└── return %12 - -julia> typeof(A[i=>1, i'=>2]) -Float64 -``` -Uh oh, that doesn't look good! Julia can't know ahead of time, based on -the inputs, what the type of the output is, besides that it will be a -`Number` (though at runtime, the output has a concrete type, `Float64`). - -So why is it designed this way? The main reason is to allow more -generic and dynamic code than traditional, statically-typed Arrays. -This allows us to have code like: -```julia -julia> i = Index(2, "i") -(dim=2|id=811|"i") - -julia> A = ITensor(i', i); - -julia> @show A; -A = ITensor ord=2 -Dim 1: (dim=2|id=811|"i")' -Dim 2: (dim=2|id=811|"i") -NDTensors.Empty{Float64,NDTensors.Dense{Float64,Array{Float64,1}}} - 2×2 - - - -julia> A[i' => 1, i => 2] = 1.2; - -julia> @show A; -A = ITensor ord=2 -Dim 1: (dim=2|id=811|"i")' -Dim 2: (dim=2|id=811|"i") -NDTensors.Dense{Float64,Array{Float64,1}} - 2×2 - 0.0 1.2 - 0.0 0.0 -``` -Here, the type of the storage of A is changed in-place. It starts as an `Empty` storage, a special trivial storage. When we set an element, we then allocate the appropriate storage. Allocations are performed only when needed, so if another element is set then no allocation is performed. -More generally, this allows ITensors to have more generic in-place -functionality, so you can write code where you don't know what the -storage is until runtime. - -This can lead to certain types of code having perfomance problems, -for example looping through ITensors with many elements can be slow: -```julia -julia> function myscale!(A::ITensor, x::Number) - for n in 1:dim(A) - A[n] = x * A[n] - end - end; - -julia> d = 10_000; - -julia> i = Index(d); - -julia> @btime myscale!(A, 2) setup = (A = random_itensor(i)); - 2.169 ms (117958 allocations: 3.48 MiB) -``` -However, this is fast: -``` -julia> function myscale!(A::Array, x::Number) - for n in 1:length(A) - A[n] = x * A[n] - end - end; - -julia> @btime myscale!(A, 2) setup = (A = randn(d)); - 3.451 μs (0 allocations: 0 bytes) - -julia> myscale2!(A::ITensor, x::Number) = myscale!(array(A), x) -myscale2! (generic function with 1 method) - -julia> @btime myscale2!(A, 2) setup = (A = random_itensor(i)); - 3.571 μs (2 allocations: 112 bytes) -``` -How does this work? It relies on a "function barrier" technique. -Julia compiles functions "just-in-time", so that calls to an inner -function written in terms of a type-stable type are still fast. -That inner function is compiled to very fast code. -The main overhead is that Julia has to determine which function -to call at runtime. - -Therefore, users should keep this in mind when they are writing -ITensors.jl code, and we warn that explicitly looping over large -ITensors by individual elements should be done with caution in -performance critical sections of your code. -However, be sure to benchmark and profile your code before -prematurely optimizing, since you may be surprised about -what are the fast and slow parts of your code. - -Some strategies for avoiding ITensor loops are: - - Use broadcasting and other built-in ITensor functionality that makes use of function barriers. - - Convert ITensors to type-stable collections like the Tensor type of NDTensors.jl and write functions in terms of the Tensor type (i.e. the function barrier techique that is used throughout ITensors.jl). - - When initializing very large ITensors elementwise, use built-in ITensor constructors, or first construct an equivalent tensor as an Array or Tensor and then convert it to an ITensor. - -## ITensor in-place operations - -In-place operations can help with optimizing code, when the -memory of the output tensor of an operation is preallocated. - -The main way to access this in ITensor is through broadcasting. -For example: -```julia -A = random_itensor(i, i') -B = random_itensor(i', i) -A .+= 2 .* B -``` -Internally, this is rewritten by Julia as a call to `broadcast!`. -ITensors.jl overloads this call (or more specifically, a lower -level function `copyto!` written in terms of a special lazy type -that saves all of the objects and operations). Then, this call is -rewritten as -```julia -map!((x,y) -> x+2*y, A, A, B) -``` -This is mostly an optimization to use when you can preallocate -storage that can be used multiple times. - -Additionally, ITensors makes the unique choice that: -```julia -C .= A .* B -``` -is interpreted as an in-place tensor contraction. What this means -is that this calls a function: -```julia -mul!(C, A, B) -``` -(likely to be given an alternative name `contract!`) which contracts -`A` and `B` into the pre-allocated memory `C`. - -Because of the design of the ITensor type (see the section above), -there is some flexibility we take in allocating memory for users. -For example, if the storage type is more narrow than the result, -for convenience we might expand it in-place. If you are worried -about memory allocations, we recommend using benchmarking and -profiling to pinpoint slow parts of your code (often times, you -may be surprised by what is actually slow). - -## NDTensors and ITensors - -ITensors.jl is built on top of another, more traditional tensor -library called NDTensors. NDTensors implements AbstractArrays with -a variety of sparse storage types, with more to come in the future. - -NDTensors implements functionality like permutation of dimensions, -fast get and set index, broadcasting, and tensor contraction (where -labels of the dimensions must be specified). - -For example: -```julia -using ITensors -using NDTensors - -T = Tensor(2,2,2) -T[1,2,1] = 1.3 # Conventional element setting - -i = Index(2) -T = Tensor((i,i',i')) # The identifiers are ignored, just interpreted as above -T[1,2,1] = 1.3 -``` -To make performant ITensor code (refer to the the previous section -on type stability and function barriers), ITensor storage data and -indices are passed by reference into Tensors, where the performance -critical operations are performed. - -An example of a function barrier using NDTensors is the following: -```julia -julia> using NDTensors - -julia> d = 10_000; - -julia> i = Index(d); - -julia> function myscale!(A::Tensor, x::Number) - for n in 1:dim(A) - A[n] = x * A[n] - end - end; - -julia> @btime myscale!(A, 2) setup = (A = Tensor(d)); - 3.530 μs (0 allocations: 0 bytes) - -julia> myscale2!(A::ITensor, x::Number) = myscale!(tensor(A), x) -myscale2! (generic function with 1 method) - -julia> @btime myscale2!(A, 2) setup = (A = random_itensor(i)); - 3.549 μs (2 allocations: 112 bytes) -``` -A very efficient function is written for the Tensor type. Then, -the ITensor version just wraps the Tensor function by calling it -after converting the ITensor to a Tensor (without any copying) -with the `tensor` function. -This is the basis for the design of all performance critical ITensors.jl functions. diff --git a/docs/src/CodeTiming.md b/docs/src/CodeTiming.md deleted file mode 100644 index 3fefa2a0a3..0000000000 --- a/docs/src/CodeTiming.md +++ /dev/null @@ -1,10 +0,0 @@ - -# Timing and Profiling your code - -It is very important to time and profile your code to make sure your code is running as fast as possible. Here are some tips on timing and profiling your code. - -If you are concerned about the performance of your code, a good place to start is Julia's [performance tips](https://docs.julialang.org/en/v1/manual/performance-tips/). - -## Timing and benchmarking - -Julia has many nice timing tools available. Tools like [@time](https://docs.julialang.org/en/v1/base/base/#Base.@time) and [TimerOutputs](https://github.com/KristofferC/TimerOutputs.jl) can be used to measure the time of specific lines of code. For microbenchmarking, we recommend the [BenchmarkTools](https://github.com/JuliaCI/BenchmarkTools.jl) package. For profiling your code, see the Julia documentation on [profiling](https://docs.julialang.org/en/v1/manual/profile/). diff --git a/docs/src/ContractionSequenceOptimization.md b/docs/src/ContractionSequenceOptimization.md index 5991946e80..183caa5e0c 100644 --- a/docs/src/ContractionSequenceOptimization.md +++ b/docs/src/ContractionSequenceOptimization.md @@ -2,7 +2,7 @@ When contracting a tensor network, the sequence of contraction makes a big difference in the computational cost. However, the complexity of determining the optimal sequence grows exponentially with the number of tensors, but there are many heuristic algorithms available for computing optimal sequences for small networks[^1][^2][^3][^4][^5][^6]. ITensors.jl provides some functionality for helping you find the optimal contraction sequence for small tensor network, as we will show below. -The algorithm in ITensors.jl currently uses a modified version of[^1] with simplifications for outer product contractions similar to those used in [TensorOperations.jl](https://github.com/Jutho/TensorOperations.jl). +The algorithm in ITensors.jl currently uses a modified version of[^1] with simplifications for outer product contractions and is based on the implementation in [TensorOperations.jl](https://github.com/Jutho/TensorOperations.jl). [^1]: [Faster identification of optimal contraction sequences for tensor networks](https://arxiv.org/abs/1304.6112) [^2]: [Improving the efficiency of variational tensor network algorithms](https://arxiv.org/abs/1310.8023) diff --git a/docs/src/DMRG.md b/docs/src/DMRG.md deleted file mode 100644 index cc1580877e..0000000000 --- a/docs/src/DMRG.md +++ /dev/null @@ -1,5 +0,0 @@ -# DMRG - -```@docs -dmrg -``` diff --git a/docs/src/DMRGObserver.md b/docs/src/DMRGObserver.md deleted file mode 100644 index 8a56cfca48..0000000000 --- a/docs/src/DMRGObserver.md +++ /dev/null @@ -1,43 +0,0 @@ -# DMRGObserver - -A DMRGObserver is a type of [observer](@ref observer) which -offers certain useful, general purpose capabilities -for DMRG calculations such as measuring custom -local observables at each step and stopping DMRG -early if certain energy convergence conditions are met. - -## Sample Usage - -In the following example, we have already made a Hamiltonian MPO `H` -and initial MPS `psi0` for a system of spins whose sites -have an associated "Sz" operator defined. We construct a -`DMRGObserver` which measures "Sz" on each site at each -step of DMRG, and also stops the calculation early if -the energy no longer changes to a relative precision of 1E-7. - -``` -Sz_observer = DMRGObserver(["Sz"],sites,energy_tol=1E-7) - -energy, psi = dmrg(H,psi0,sweeps,observer=Sz_observer) - -for (sw,Szs) in enumerate(measurements(Sz_observer)["Sz"]) - println("Total Sz after sweep $sw = ", sum(Szs)/N) -end -``` - - -## Constructors - -```@docs -DMRGObserver(;energy_tol::Float64,minsweeps::Int) -DMRGObserver(ops::Vector{String},sites::Vector{<:Index};energy_tol::Float64,minsweeps::Int) -``` - -## Methods - -```@docs -measurements(::DMRGObserver) -DMRGMeasurement -energies(::DMRGObserver) -``` - diff --git a/docs/src/DeveloperGuide.md b/docs/src/DeveloperGuide.md deleted file mode 100644 index c30269e69e..0000000000 --- a/docs/src/DeveloperGuide.md +++ /dev/null @@ -1,100 +0,0 @@ -# Developer Guide - - -## Keyword Argument Best Practices - -Keyword arguments such as `f(x,y; a=1, b=2)` are a powerful Julia feature, but it is easy to -misuse them in library code. Below are the "best practices" for using keyword arguments -when developing ITensor library code. - -A particular challenge how to properly use keyword argument "forwarding" where the notation -`f(; a, b, kwargs...)` allows any number of keyword arguments to be passed. -If a keyword argument is misspelled, then forwarding keywords with `kwargs...` will -silently allow the misspelling, whereas ideally there would be an error message. - -Best practices: - -1. **Popping Terminal Keyword Arguments**: - When passing keyword arguments downward through a stack of function calls, if a certain keyword - argument will not be used in any functions further down the stack, then these arguments should - be *listed explicitly* to remove them from the keyword arguments. - - For example, in a call stack `fA -> fB -> fC` if a keyword argument - such as `cutoff` is used in the body of `fB` but not in `fC`, then use the following pattern: - - ``` - function fA(...; kwargs...) - ... - fB(...; kwargs...) - ... - end - - function fB(...; cutoff, kwargs...) # <- explicitly list cutoff here - ... - truncate!(psi; cutoff) # <- fB uses cutoff - fC(...; kwargs...) # fC does not get passed cutoff - end - - function fC(...; maxdim, outputlevel) # fC does not use or need the `cutoff` kwarg - ... - end - ``` - -2. **Leaf Functions Should Not Take `kwargs...`**: - Functions which are the last in the call stack to take any keyword arguments - should not take keyword arguments by the `kwargs...` pattern. They should only take an explicit - list of keyword arguments, so as to ensure that an error is thrown if a keyword argument - is misspelled or missing (if it has no default value). - - Example: `fC` above is a leaf function and does not have `kwargs...` in its signature. - -3. **Use Functions to Set Defaults**: - Keyword arguments can be made optional by providing default values. To avoid having explicit and - possibly inconsistent defaults spread all over the library code, use globally defined functions to - provide these defaults. - - For example: - ``` - function sum(A::MPS, B::MPS; cutoff=default_cutoff(), kwargs...) - ... - end - - function inner(A::MPS, B::MPS; cutoff=default_cutoff(), kwargs...) - ... - end - ``` - where above the default value for the `cutoff` keyword is provided by a function `default_cutoff()` - that is defined for the whole library. - -4. **Use Named Tuples to "Tunnel" Keywords to Leaf Functions**: - This is a more advanced pattern. In certain situations, there might be multiple leaf - functions depending on the execution pathway of the code or in cases where the leaf function - is a "callback" passed into the code from the upper-level calling code. - - In such cases, different leaf function implementations may expect different sets of keyword arguments. - - To avoid requiring all leaf functions to take all possible keyword arguments (or to use the `kwargs...` - pattern as a workaround, breaking rule #2 above), use the following pattern: - - ``` - function fA(callback, psi; callback_args, kwargs...) - ... - callback(psi; callback_args...) - ... - end - - my_callback(psi; a, b) = ... # define custom callback function - - # Call fA like this: - fA(my_callback, psi; callback_args = (; a, b)) - - ``` - -5. **External (non-ITensor) Functions**: - Though it requires judgment in each case, if the keyword arguments an external - (non-ITensor) function accepts are small in number, not expected to change, - and known ahead of time, try to list them explicitly if possible (rather than forwarding - with `kwargs...`). Possible exceptions could be if you want to make use of defaults - defined for keyword arguments of an external function. - - diff --git a/docs/src/Einsum.md b/docs/src/Einsum.md deleted file mode 100644 index 873ac3239d..0000000000 --- a/docs/src/Einsum.md +++ /dev/null @@ -1,192 +0,0 @@ -# ITensor Index identity: dimension labels and Einstein notation - -Many tensor contraction libraries use [Einstein notation](https://en.wikipedia.org/wiki/Einstein_notation), -such as [NumPy's einsum function](https://numpy.org/doc/stable/reference/generated/numpy.einsum.html), [ncon](https://arxiv.org/abs/1402.0939), and various Julia packages such as [TensorOperations.jl](https://github.com/Jutho/TensorOperations.jl), [Tullio.jl](https://github.com/mcabbott/Tullio.jl), [OMEinsum.jl](https://github.com/under-Peter/OMEinsum.jl), and [Einsum.jl](https://github.com/ahwillia/Einsum.jl), among others. - -ITensor also uses Einstein notation, however the labels are stored inside the tensor and carried around with them during various operations. In addition, the labels that determine if tensor indices match with each other, and therefore automatically contract when doing `*` or match when adding or subtracting, are more sophisticated than simple characters or strings. ITensor indices are given a unique random ID number when they are constructed, and additionally users can add additional information like prime levels and tags which uniquely determine an Index. This is in contrast to simpler implementations of the same idea, such as the [NamedDims.jl](https://github.com/invenia/NamedDims.jl) package, which only allow symbols as the metadata for uniquely identifying a tensor/array dimension. - -```@setup itensor -using ITensors -using Random -Random.seed!(1) -``` - -## Index identity - -Here is an illustration of how the different types of Index metadata (random ID, prime level, and tags) work for Index identity: -```@repl itensor -i = Index(2) -j = Index(2) -i == j -id(i) -id(j) -ip = i' -ip == i -plev(i) == 0 -plev(ip) == 1 -noprime(ip) == i -ix = addtags(i, "x") -ix == i -removetags(ix, "x") == i -ixyz = addtags(ix, "y,z") -ixyz == addtags(i, "z,y,x") -``` - -The different metadata that are stored inside of ITensor indices that determine their identity are useful in different contexts. The random ID is particularly useful in the case when a new Index needs to be generated internally by ITensor, such as when performing a matrix factorization. In the case of a matrix factorization, we want to make sure that the new Index will not accidentally clash with an existing one, for example: -```@repl itensor -i = Index(2, "i") -j = Index(2, "j") -A = random_itensor(i, j) -U, S, V = svd(A, i; lefttags="i", righttags="j"); -inds(U) -inds(S) -inds(V) -norm(U * S * V - A) -``` -You can see that it would have been a problem here if there wasn't a new ID assigned to the Index, since it would have clashed with the original index. In this case, it could be avoided by giving the new indices different tags (with the keyword arguments `lefttags` and `righttags`), but in more complicated examples where it is not practical to do that (such as a case where many new indices are being introduced, for example for a tensor train (TT)/matrix product state (MPS)), it is convenient to not force users to come up with unique prime levels or tags themselves. It can also help to avoid accidental contractions in more complicated tensor network algorithms where there are many indices that can potentially have the same prime levels or tags. - -In contrast, using multiple indices with the same Index ID but different prime levels and tags can be useful in situations where there is a more fundamental relationship between the spaces. For example, in the case of an ITensor corresponding to a Hermitian operator, it is helpful to make the bra space and ket spaces the same up to a prime level: -```repl itensor -i = Index(2, "i") -j = Index(3, "j") -A = random_itensor(i', j', dag(i), dag(j)) -H = 0.5 * (A + swapprime(dag(A), 0 => 1)) -v = random_itensor(i, j) -Hv = noprime(H * v) -vH = dag(v)' * H -norm(Hv - dag(vH)) -``` -Note that we have added `dag` in a few places, which is superfluous in this case since the tensors are real and dense but become important when the tensors are complex and/or have symmetries. -You can see that in this case, it is very useful to relate the bra and ket spaces by prime levels, since it makes it much easier to perform operations that map from one space to another. We could have created `A` from 4 entirely different indices with different ID numbers, but it would make the operations a bit more cumbersome, as shown below: -```@repl itensor -i = Index(2, "i") -j = Index(3, "j") -ip = Index(2, "i") -jp = Index(3, "jp") -A = random_itensor(ip, jp, dag(i), dag(j)) -H = 0.5 * (A + swapinds(dag(A), (i, j), (ip, jp))) -v = random_itensor(i, j) -Hv = replaceinds(H * v, (ip, jp) => (i, j)) -vH = replaceinds(dag(v), (i, j) => (ip, jp)) * H -norm(Hv - dag(vH)) -``` - -## Relationship to other Einstein notation-based libraries - -Here we show examples of different ways to perform the contraction -`"ab,bc,cd->ad"` in ITensor. - -```@repl itensor -da, dc = 2, 3; -db, dd = da, dc; -tags = ("a", "b", "c", "d"); -dims = (da, db, dc, dd); -a, b, c, d = Index.(dims, tags); -Aab = random_itensor(a, b) -Bbc = random_itensor(b, c) -Ccd = random_itensor(c, d) - -# "ab,bc,cd->ad" -out1 = Aab * Bbc * Ccd -@show hassameinds(out1, (a, d)) - -# -# Using replaceinds (most general way) -# - -# "ba,bc,dc->ad" -Aba = replaceinds(Aab, (a, b) => (b, a)) -Cdc = replaceinds(Ccd, (c, d) => (d, c)) -out2 = Aba * Bbc * Cdc -@show hassameinds(out2, (a, d)) - -# -# Using setinds -# - -# This is a bit lower level -# since it doesn't check if the indices -# are compatible in dimension, -# so is not recommended in general. -using ITensors: setinds - -Aba = setinds(Aab, (b, a)) -Cdc = setinds(Ccd, (d, c)) -out2 = Aba * Bbc * Cdc -@show hassameinds(out2, (a, d)) - -# -# Using prime levels (assuming -# the indices were made with these -# prime levels in the first place) -# - -a = Index(da, "a") -c = Index(dc, "c") -b, d = a', c' -Aab = random_itensor(a, b) -Bbc = random_itensor(b, c) -Ccd = random_itensor(c, d) -out1 = Aab * Bbc * Ccd -@show hassameinds(out1, (a, d)) - -Aba = swapprime(Aab, 0 => 1) -Cdc = swapprime(Ccd, 0 => 1) -out2 = Aba * Bbc * Cdc -@show hassameinds(out2, (a, d)) - -# -# Using tags (assuming -# the indices were made with these -# tags in the first place) -# - -a = Index(da, "a") -c = Index(dc, "c") -b, d = settags(a, "b"), settags(c, "d") -Aab = random_itensor(a, b) -Bbc = random_itensor(b, c) -Ccd = random_itensor(c, d) -out1 = Aab * Bbc * Ccd -@show hassameinds(out1, (a, d)) - -Aba = swaptags(Aab, "a", "b") -Cdc = swaptags(Ccd, "c", "d") -out2 = Aba * Bbc * Cdc -@show hassameinds(out2, (a, d)) - -# -# Using Julia Arrays -# - -A = randn(da, db) -B = randn(db, dc) -C = randn(dc, dd) - -tags = ("a", "b", "c", "d") -dims = (da, db, dc, dd) -a, b, c, d = Index.(dims, tags) - -Aab = itensor(A, a, b) -Bbc = itensor(B, b, c) -Ccd = itensor(C, c, d) -out1 = Aab * Bbc * Ccd -@show hassameinds(out1, (a, d)) - -Aba = itensor(A, b, a) -Cdc = itensor(C, d, c) -out2 = Aba * Bbc * Cdc -@show hassameinds(out2, (a, d)) - -# -# Note that we may start allowing -# this notation in future: -# (https://github.com/ITensor/ITensors.jl/issues/673) -# -#out1 = A[a, b] * B[b, c] * C[c, d] -#@show hassameinds(out1, (a, d)) -# -#out2 = A[b, a] * B[b, c] * C[d, c] -#@show hassameinds(out2, (a, d)) -``` - diff --git a/docs/src/HDF5FileFormats.md b/docs/src/HDF5FileFormats.md index 47cf26ea8f..855c730070 100644 --- a/docs/src/HDF5FileFormats.md +++ b/docs/src/HDF5FileFormats.md @@ -154,34 +154,3 @@ Datasets and Subgroups: by the HDF5.jl library for storing `Vector{Float64}` or `Vector{ComplexF64}`) -## [MPS](@id mps_hdf5) - -HDF5 file format for `ITensorMPS.MPS` - -Attributes: -* "version" = 1 -* "type" = "MPS" - -Datasets and Subgroups: -* "length" [integer] = number of tensors of the MPS -* "rlim" [integer] = right orthogonality limit -* "llim" [integer] = left orthogonality limit -* "MPS[n]" [group,ITensor] = each of these groups, where n=1,...,length, stores the nth ITensor of the MPS - - -## [MPO](@id mpo_hdf5) - -HDF5 file format for `ITensorMPS.MPO` - -Attributes: -* "version" = 1 -* "type" = "MPO" - -Datasets and Subgroups: -* "length" [integer] = number of tensors of the MPO -* "rlim" [integer] = right orthogonality limit -* "llim" [integer] = left orthogonality limit -* "MPO[n]" [group,ITensor] = each of these groups, where n=1,...,length, stores the nth ITensor of the MPO - - - diff --git a/docs/src/IncludedSiteTypes.md b/docs/src/IncludedSiteTypes.md deleted file mode 100644 index f8396c8416..0000000000 --- a/docs/src/IncludedSiteTypes.md +++ /dev/null @@ -1,388 +0,0 @@ -# SiteTypes Included with ITensor - -## "S=1/2" SiteType - -Site indices with the "S=1/2" site type represent ``S=1/2`` spins with the states -``|\!\uparrow\rangle``, ``|\!\downarrow\rangle``. - -Making a single "S=1/2" site or collection of N "S=1/2" sites -``` -s = siteind("S=1/2") -sites = siteinds("S=1/2",N) -``` - -Available keyword arguments for enabling and customizing quantum numbers (QN) subspaces: -- `conserve_qns` (default: false): conserve total ``S^z`` -- `conserve_sz` (default: conserve_qns): conserve total ``S^z`` -- `conserve_szparity` (default: false): conserve total ``S^z`` modulo two -- `qnname_sz` (default: "Sz"): name of total ``S^z`` QN -- `qnname_szparity` (default: "SzParity"): name of total ``S^z`` modulo two QN -For example: -``` -sites = siteinds("S=1/2",N; conserve_szparity=true, qnname_szparity="SzP") -``` - -Operators associated with "S=1/2" sites can be made using the `op` function, -for example -``` -Sz = op("Sz",s) -Sz4 = op("Sz",sites[4]) -``` - -Available operators are exactly the same as those for the "Qubit" site type. Please -see the list of "Qubit" operators below. - -## "Qubit" SiteType - -Site indices with the "Qubit" site type represent qubits with the states -``|0\rangle``, ``|1\rangle``. - -Making a single "Qubit" site or collection of N "Qubit" sites -``` -s = siteind("Qubit") -sites = siteinds("Qubit",N) -``` - -Available keyword arguments for enabling and customizing quantum numbers (QN) subspaces: -- `conserve_qns` (default: false): conserve total qubit parity -- `conserve_parity` (default: conserve_qns): conserve total qubit parity -- `conserve_number` (default: false): conserve total qubit number -- `qnname_parity` (default: "Parity"): name of total qubit parity QN -- `qnname_number` (default: "Number"): name of total qubit number QN -For example: -``` -sites = siteinds("Qubit",N; conserve_parity=true) -``` - -#### "Qubit" and "S=1/2" States - -The available state names for "Qubit" sites are: -- `"0"` (aliases: `"Z+"`, `"Zp"`, `"Up"`, `"↑"`) Qubit in the 0 state -- `"1"` (aliases: `"Z-"`, `"Zm"`, `"Dn"`, `"↓"`) Qubit in the 1 state -- `"+"` (aliases: `"X+"`, `"Xp"`) Qubit in the $|+\rangle$ state (+1 eigenvector of $\sigma_x$) -- `"+"` (aliases: `"X-"`, `"Xm"`) Qubit in the $|-\rangle$ state (-1 eigenvector of $\sigma_x$) -- `"i"` (aliases: `"Y+"`, `"Yp"`) Qubit in the $|i\rangle$ state (+1 eigenvector of $\sigma_y$) -- `"-i"` (aliases: `"Y-"`, `"Ym"`) Qubit in the $|-i\rangle$ state (+1 eigenvector of $\sigma_y$) - -#### "Qubit" and "S=1/2" Operators - -Operators or gates associated with "Qubit" sites can be made using the `op` function, -for example -``` -H = op("H",s) -H3 = op("H",sites[3]) -``` - -Single-qubit operators: -- `"X"` (aliases: `"σx"`, `"σ1"`) Pauli X operator -- `"Y"` (aliases: `"σy"`, `"σ2"`) Pauli Y operator -- `"iY"` (aliases: `"iσy"`, `"iσ2"`) Pauli Y operator times i -- `"Z"` (aliases: `"σz"`, `"σ3"`) Pauli Z operator -- `"√NOT"` (aliases: `"X"`) -- `"H"` Hadamard gate -- `"Phase"` (takes optional argument: ϕ=π/2) (aliases: `"P"`, `"S"`) -- `"π/8"` (aliases: `"T"`) -- `"Rx"` (takes argument: θ) Rotation around x axis -- `"Ry"` (takes argument: θ) Rotation around y axis -- `"Rz"` (takes argument: θ) Rotation around z axis -- `"Rn"` (takes arguments: θ, ϕ, λ) (aliases: `"Rn̂"`) Rotation about axis n=(θ, ϕ, λ) -- `"Proj0"` (aliases: `"ProjUp"`, `"projUp"`) Operator $|0\rangle\langle 0|$ -- `"Proj1"` (aliases: `"ProjDn"`, `"projDn"`) Operator $|1\rangle\langle 1|$ - -Spin operators: -- `"Sz"` (aliases: `"Sᶻ"`) Spin z operator $S^z = \frac{1}{2} \sigma_z$ -- `"S+"` (alises: `"S⁺"`, `"Splus"`) Raising operator $S^+ = S^x + iS^y$ -- `"S-"` (aliases: `"S⁻"`, `"Sminus"`) Lowering operator $S^- = S^x - iS^y$ -- `"Sx"` (alises: `"Sˣ"`) Spin x operator $S^x = \frac{1}{2} \sigma_x$ -- `"iSy"` (aliases: `"iSʸ"`) i times spin y operator $iS^y = \frac{i}{2} \sigma_y$ -- `"Sy"` (aliases: `"Sʸ"`) Spin y operator $S^y = \frac{1}{2} \sigma_y$ -- `"S2"` (aliases: "S²"`) Square of spin vector operator $S^2=\vec{S}\cdot\vec{S}=\frac{3}{4} I$ -- `"ProjUp"` (aliases: `"projUp"`, `"Proj0"`) Operator $|\!↑\rangle\langle ↑\!|$ -- `"ProjDn"` (aliases: `"projDn"`, `"Proj1"`) Operator $|\!↓\rangle\langle ↓\!|$ - -Two-qubit gates: -- `"CNOT"` (aliases: `"CX"`) Controlled NOT gate -- `"CY"` Controlled Y gate -- `"CZ"` Controlled Z gate -- `"CPHASE"` (aliases: `"Cphase"`) Controlled Phase gate -- `"CRx"` (aliases: `"CRX"`) (takes arguments: θ) -- `"CRy"` (aliases: `"CRY"`) (takes arguments: θ) -- `"CRz"` (aliases: `"CRZ"`) (takes arguments: θ) -- `"CRn"` (aliases: `"CRn̂"`) (takes arguments: θ, ϕ, λ) -- `"SWAP"` (aliases: `"Swap"`) -- `"√SWAP"` (aliases: `"√Swap"`) -- `"iSWAP"` (aliases: `"iSwap"`) -- `"√iSWAP"` (aliases: `"√iSwap"`) -- `"Rxx"` (aliases: `"RXX"`) (takes arguments: ϕ) Ising (XX) coupling gate -- `"Ryy"` (aliases: `"RYY"`) (takes arguments: ϕ) Ising (YY) coupling gate -- `"Rzz"` (aliases: `"RZZ"`) (takes arguments: ϕ) Ising (ZZ) coupling gate - -Three-qubit gates: -- `"Toffoli"` (aliases `"CCNOT"`, `"CCX"`, `"TOFF"`) -- `"Fredkin"` (aliases `"CSWAP"`, `"CSwap"`, `"CS"`) - -Four-qubit gates: -- `"CCCNOT"` - -## "S=1" SiteType - -Site indices with the "S=1" site type represent ``S=1`` spins with the states -``|\!\uparrow\rangle``, ``|0\rangle``, ``|\!\downarrow\rangle``. - -Making a single "S=1" site or collection of N "S=1" sites -``` -s = siteind("S=1") -sites = siteinds("S=1",N) -``` - -Available keyword arguments for enabling and customizing quantum numbers (QN) subspaces: -- `conserve_qns` (default: false): conserve total ``S^z`` -- `conserve_sz` (default: conserve_qns): conserve total ``S^z`` -- `qnname_sz` (default: "Sz"): name of total ``S^z`` QN -For example: -``` -sites = siteinds("S=1",N; conserve_sz=true, qnname_sz="TotalSz") -``` - -#### "S=1" States - -The available state names for "S=1" sites are: -- `"Up"` (aliases: `"Z+"`, `"↑"`) spin in the up state -- `"Z0"` (aliases: `"0"`) spin in the Sz=0 state -- `"Dn"` (aliases: `"Z-"`, `"↓"`) spin in the down state - -#### "S=1" Operators - -Operators associated with "S=1" sites can be made using the `op` function, -for example -``` -Sz = op("Sz",s) -Sz4 = op("Sz",sites[4]) -``` - -Spin operators: -- `"Sz"` (aliases: `"Sᶻ"`) -- `"Sz2"` Square of `S^z` operator -- `"S+"` (alises: `"S⁺"`, `"Splus"`) -- `"S-"` (aliases: `"S⁻"`, `"Sminus"`) -- `"Sx"` (alises: `"Sˣ"`) -- `"Sx2"` Square of `S^x` operator -- `"iSy"` (aliases: `"iSʸ"`) -- `"Sy"` (aliases: `"Sʸ"`) -- `"Sy2"` Square of `S^y` operator -- `"S2"` (aliases: "S²"`) - -## "Boson" SiteType - -The "Boson" site type is an alias for the "Qudit" site type. Please -see more information about "Qudit" below: - -## "Qudit" SiteType - -Making a single "Qudit" site or collection of N "Qudit" sites -``` -s = siteind("Qudit") -sites = siteinds("Qudit",N) -``` - -Available keyword arguments for enabling and customizing quantum numbers (QN) subspaces: -- `dim` (default: 2): dimension of the index (number of qudit or boson values) -- `conserve_qns` (default: false): conserve total qudit or boson number -- `conserve_number` (default: conserve_qns): conserve total qudit or boson number -- `qnname_number` (default: "Number"): name of total qudit or boson number QN -For example: -``` -sites = siteinds("Qudit",N; conserve_number=true) -``` - -#### "Qudit" and "Boson" Operators - -Operators associated with "Qudit" sites can be made using the `op` function, -for example -``` -A = op("A",s) -A4 = op("A",sites[4]) -``` - -Single-qudit operators: -- `"A"` (aliases: `"a"`) -- `"Adag"` (aliases: `"adag"`, `"a†"`) -- `"N"` (aliases: `"n"`) - -Two-qudit operators: -- `"ab"` -- `"a†b"` -- `"ab†"` -- `"a†b†"` - -## "Fermion" SiteType - -Site indices with the "Fermion" SiteType represent -spinless fermion sites with the states -``|0\rangle``, ``|1\rangle``, corresponding to zero fermions or one fermion. - -Making a single "Fermion" site or collection of N "Fermion" sites -``` -s = siteind("Fermion") -sites = siteinds("Fermion",N) -``` - -Available keyword arguments for enabling and customizing quantum numbers (QN) subspaces: -- `conserve_qns` (default: false): conserve total number of fermions -- `conserve_nf` (default: conserve_qns): conserve total number of fermions -- `conserve_nfparity` (default: conserve_qns): conserve total fermion number parity -- `qnname_nf` (default: "Nf"): name of total fermion number QN -- `qnname_nfparity` (default: "NfParity"): name of total fermion number parity QN -For example: -``` -sites = siteinds("Fermion",N; conserve_nfparity=true) -``` - -#### "Fermion" States - -The available state names for "Fermion" sites are: -- `"0"` (aliases: `"Emp"`) unoccupied fermion site -- `"1"` (aliases: `"Occ"`) occupied fermion site - -#### "Fermion" Operators - -Operators associated with "Fermion" sites can be made using the `op` function, -for example -``` -C = op("C",s) -C4 = op("C",sites[4]) -``` - -Single-fermion operators: -- `"N"` (aliases: `"n"`) Density operator -- `"C"` (aliases: `"c"`) Fermion annihilation operator -- `"Cdag"` (aliases: `"cdag"`, `"c†"`) Fermion creation operator -- `"F"` Jordan-Wigner string operator - -## "Electron" SiteType - -The states of site indices with the "Electron" SiteType correspond to -``|0\rangle``, ``|\!\uparrow\rangle``, ``|\!\downarrow\rangle``, ``|\!\uparrow\downarrow\rangle``. - -Making a single "Electron" site or collection of N "Electron" sites -``` -s = siteind("Electron") -sites = siteinds("Electron",N) -``` - -Available keyword arguments for enabling and customizing quantum numbers (QN) subspaces: -- `conserve_qns` (default: false): conserve total number of electrons -- `conserve_sz` (default: conserve_qns): conserve total ``S^z`` -- `conserve_nf` (default: conserve_qns): conserve total number of electrons -- `conserve_nfparity` (default: conserve_qns): conserve total electron number parity -- `qnname_sz` (default: "Sz"): name of total ``S^z`` QN -- `qnname_nf` (default: "Nf"): name of total electron number QN -- `qnname_nfparity` (default: "NfParity"): name of total electron number parity QN -For example: -``` -sites = siteinds("Electron",N; conserve_nfparity=true) -``` - -#### "Electron" States - -The available state names for "Electron" sites are: -- `"Emp"` (aliases: `"0"`) unoccupied electron site -- `"Up"` (aliases: `"↑"`) electron site occupied with one up electron -- `"Dn"` (aliases: `"↓"`) electron site occupied with one down electron -- `"UpDn"` (aliases: `"↑↓"`) electron site occupied with two electrons (one up, one down) - -#### "Electron" Operators - -Operators associated with "Electron" sites can be made using the `op` function, -for example -``` -Cup = op("Cup",s) -Cup4 = op("Cup",sites[4]) -``` - -Single-fermion operators: -- `"Ntot"` (aliases: `"ntot"`) Total density operator -- `"Nup"` (aliases: `"n↑"`) Up density operator -- `"Ndn"` (aliases: `"n↓"`) Down density operator -- `"Cup"` (aliases: `"c↑"`) Up-spin annihilation operator -- `"Cdn"` (aliases: `"c↓"`) Down-spin annihilation operator -- `"Cdagup"` (aliases: `"c†↑"`) Up-spin creation operator -- `"Cdagdn"` (aliases: `"c†↓"`) Down-spin creation operator -- `"Sz"` (aliases: `"Sᶻ"`) -- `"Sx"` (aliases: `"Sˣ"`) -- `"S+"` (aliases: `"Sp"`, `"S⁺"`,`"Splus"`) -- `"S-"` (aliases: `"Sm"`, `"S⁻"`, `"Sminus"`) -- `"F"` Jordan-Wigner string operator -- `"Fup"` (aliases: `"F↑"`) Up-spin Jordan-Wigner string operator -- `"Fdn"` (aliases: `"F↓"`) Down-spin Jordan-Wigner string operator - -Non-fermionic single particle operators (these do not have Jordan-Wigner string attached, -so will commute within systems such as OpSum or the `apply` function): -- `"Aup"` (aliases: `"a↑"`) Up-spin annihilation operator -- `"Adn"` (aliases: `"a↓"`) Down-spin annihilation operator -- `"Adagup"` (aliases: `"a†↑"`) Up-spin creation operator -- `"Adagdn"` (aliases: `"a†↓"`) Down-spin creation operator - - -## "tJ" SiteType - -"tJ" sites are similar to electron sites, but cannot be doubly occupied -The states of site indices with the "tJ" SiteType correspond to -``|0\rangle``, ``|\!\uparrow\rangle``, ``|\!\downarrow\rangle``. - -Making a single "tJ" site or collection of N "tJ" sites -``` -s = siteind("tJ") -sites = siteinds("tJ",N) -``` - -Available keyword arguments for enabling and customizing quantum numbers (QN) subspaces: -- `conserve_qns` (default: false): conserve total number of fermions -- `conserve_nf` (default: conserve_qns): conserve total number of fermions -- `conserve_nfparity` (default: conserve_qns): conserve total fermion number parity -- `qnname_nf` (default: "Nf"): name of total fermion number QN -- `qnname_nfparity` (default: "NfParity"): name of total fermion number parity QN -For example: -``` -sites = siteinds("tJ",N; conserve_nfparity=true) -``` - -#### "tJ" States - -The available state names for "tJ" sites are: -- `"Emp"` (aliases: `"0"`) unoccupied site -- `"Up"` (aliases: `"↑"`) site occupied with one up electron -- `"Dn"` (aliases: `"↓"`) site occupied with one down electron - -#### "tJ" Operators - -Operators associated with "tJ" sites can be made using the `op` function, -for example -``` -Cup = op("Cup",s) -Cup4 = op("Cup",sites[4]) -``` - -Single-fermion operators: -- `"Ntot"` (aliases: `"ntot"`) Total density operator -- `"Nup"` (aliases: `"n↑"`) Up density operator -- `"Ndn"` (aliases: `"n↓"`) Down density operator -- `"Cup"` (aliases: `"c↑"`) Up-spin annihilation operator -- `"Cdn"` (aliases: `"c↓"`) Down-spin annihilation operator -- `"Cdagup"` (aliases: `"c†↑"`) Up-spin creation operator -- `"Cdagdn"` (aliases: `"c†↓"`) Down-spin creation operator -- `"Sz"` (aliases: `"Sᶻ"`) -- `"Sx"` (aliases: `"Sˣ"`) -- `"S+"` (aliases: `"Sp"`, `"S⁺"`,`"Splus"`) -- `"S-"` (aliases: `"Sm"`, `"S⁻"`, `"Sminus"`) -- `"F"` Jordan-Wigner string operator -- `"Fup"` (aliases: `"F↑"`) Up-spin Jordan-Wigner string operator -- `"Fdn"` (aliases: `"F↓"`) Down-spin Jordan-Wigner string operator - -Non-fermionic single particle operators (these do not have Jordan-Wigner string attached, -so will commute within systems such as OpSum or the `apply` function): -- `"Aup"` (aliases: `"a↑"`) Up-spin annihilation operator -- `"Adn"` (aliases: `"a↓"`) Down-spin annihilation operator -- `"Adagup"` (aliases: `"a†↑"`) Up-spin creation operator -- `"Adagdn"` (aliases: `"a†↓"`) Down-spin creation operator - diff --git a/docs/src/IndexSetType.md b/docs/src/IndexSetType.md deleted file mode 100644 index ce69a81d2c..0000000000 --- a/docs/src/IndexSetType.md +++ /dev/null @@ -1,32 +0,0 @@ -# Index collections - -Collections of `Index` are used throughout ITensors.jl to represent the dimensions of tensors. In general, collections that are recognized and returned by ITensors.jl functions are either `Vector` of `Index` or `Tuple` of `Index`, depending on the context. For example internally an `ITensor` has a static number of indices so stores a `Tuple` of `Index`, while set operations like `commoninds((i, j, k), (j, k, l))` will return a `Vector` `[j, k]` since the operation is inherently dynamic, i.e. the number of indices in the intersection can't in general be known before running the code. `Vector` of `Index` and `Tuple` of `Index` can usually be used interchangeably, but one or the other may be faster depending on the operation being performed. - -## [Priming and tagging](@id Priming_and_tagging_IndexSet) - -Documentation for priming and tagging collections of Index can be found in the ITensor [Priming and tagging](@ref Priming_and_tagging_ITensor) section. - -## Set operations - -Documentation for set operations involving Index collections can be found in the ITensor [Index collections set operations](@ref) section. - -## Subsets - -```@docs -getfirst(::Function, ::IndexSet) -getfirst(::IndexSet) -``` - -## Iterating - -```@docs -eachval(::Index...) -eachindval(::Index...) -``` - - -## Symmetry related properties - -```@docs -dir(::IndexSet, ::Index) -``` diff --git a/docs/src/MPSandMPO.md b/docs/src/MPSandMPO.md deleted file mode 100644 index 5474bbc1ae..0000000000 --- a/docs/src/MPSandMPO.md +++ /dev/null @@ -1,162 +0,0 @@ -# MPS and MPO - -## Types - -```@docs -MPS -MPO -``` - -## MPS Constructors - -```@docs -MPS(::Int) -MPS(::Type{<:Number}, ::Vector{<:Index}) -random_mps(sites::Vector{<:Index}) -random_mps(::Type{<:Number}, sites::Vector{<:Index}) -random_mps(::Vector{<:Index}, ::Any) -MPS(::Vector{<:Index}, ::Any) -MPS(::Type{<:Number}, ::Vector{<:Index}, ::Any) -MPS(::Vector{<:Pair{<:Index}}) -MPS(::Type{<:Number}, ::Vector{<:Pair{<:Index}}) -``` - -## MPO Constructors - -```@docs -MPO(::Int) -MPO(::Type{<:Number}, ::Vector{<:Index}, ::Vector{String}) -MPO(::Type{<:Number}, ::Vector{<:Index}, ::String) -``` - -## Copying behavior - -```@docs -copy(::ITensorMPS.AbstractMPS) -deepcopy(::ITensorMPS.AbstractMPS) -``` - -## Properties - -```@docs -eltype(::ITensorMPS.AbstractMPS) -flux(::ITensorMPS.AbstractMPS) -hasqns(::ITensorMPS.AbstractMPS) -length(::ITensorMPS.AbstractMPS) -maxlinkdim(::ITensorMPS.AbstractMPS) -``` - -## Obtaining and finding indices - -```@docs -siteinds(::typeof(commoninds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS, ::Int) -siteinds(::typeof(uniqueinds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS, ::Int) -findsite -findsites -firstsiteinds -linkind(::ITensorMPS.AbstractMPS,::Int) -siteind(::MPS, ::Int) -siteind(::typeof(first), ::MPS, ::Int) -siteinds(::MPS) -siteind(::MPO, ::Int) -siteinds(::MPO) -siteinds(::ITensorMPS.AbstractMPS, ::Int) -``` - -## Priming and tagging - -```@docs -prime(::ITensorMPS.AbstractMPS) -prime(::typeof(siteinds), ::ITensorMPS.AbstractMPS) -prime(::typeof(linkinds), ::ITensorMPS.AbstractMPS) -prime(::typeof(siteinds), ::typeof(commoninds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) -prime(::typeof(siteinds), ::typeof(uniqueinds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) - -swapprime(::ITensorMPS.AbstractMPS, args...; kwargs...) - -setprime(::ITensorMPS.AbstractMPS) -setprime(::typeof(siteinds), ::ITensorMPS.AbstractMPS) -setprime(::typeof(linkinds), ::ITensorMPS.AbstractMPS) -setprime(::typeof(siteinds), ::typeof(commoninds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) -setprime(::typeof(siteinds), ::typeof(uniqueinds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) - -noprime(::ITensorMPS.AbstractMPS) -noprime(::typeof(siteinds), ::ITensorMPS.AbstractMPS) -noprime(::typeof(linkinds), ::ITensorMPS.AbstractMPS) -noprime(::typeof(siteinds), ::typeof(commoninds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) -noprime(::typeof(siteinds), ::typeof(uniqueinds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) - -addtags(::ITensorMPS.AbstractMPS) -addtags(::typeof(siteinds), ::ITensorMPS.AbstractMPS) -addtags(::typeof(linkinds), ::ITensorMPS.AbstractMPS) -addtags(::typeof(siteinds), ::typeof(commoninds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) -addtags(::typeof(siteinds), ::typeof(uniqueinds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) - -removetags(::ITensorMPS.AbstractMPS) -removetags(::typeof(siteinds), ::ITensorMPS.AbstractMPS) -removetags(::typeof(linkinds), ::ITensorMPS.AbstractMPS) -removetags(::typeof(siteinds), ::typeof(commoninds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) -removetags(::typeof(siteinds), ::typeof(uniqueinds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) - -replacetags(::ITensorMPS.AbstractMPS) -replacetags(::typeof(siteinds), ::ITensorMPS.AbstractMPS) -replacetags(::typeof(linkinds), ::ITensorMPS.AbstractMPS) -replacetags(::typeof(siteinds), ::typeof(commoninds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) -replacetags(::typeof(siteinds), ::typeof(uniqueinds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) - -settags(::ITensorMPS.AbstractMPS) -settags(::typeof(siteinds), ::ITensorMPS.AbstractMPS) -settags(::typeof(linkinds), ::ITensorMPS.AbstractMPS) -settags(::typeof(siteinds), ::typeof(commoninds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) -settags(::typeof(siteinds), ::typeof(uniqueinds), ::ITensorMPS.AbstractMPS, ::ITensorMPS.AbstractMPS) -``` - -## Operations - -```@docs -expect(::MPS, ::Any) -correlation_matrix(::MPS, ::AbstractString, ::AbstractString) -dag(::ITensorMPS.AbstractMPS) -dense(::ITensorMPS.AbstractMPS) -movesite(::ITensorMPS.AbstractMPS, ::Pair{Int, Int};orthocenter::Int,kwargs...) -orthogonalize! -replacebond!(::MPS, ::Int, ::ITensor) -sample(::MPS) -sample!(::MPS) -sample(::MPO) -swapbondsites(::ITensorMPS.AbstractMPS, ::Int; kwargs...) -truncate! -``` - -## Gate evolution - -```@docs -product(::ITensor, ::ITensorMPS.AbstractMPS) -product(::Vector{ITensor}, ::ITensorMPS.AbstractMPS) -``` - -## Algebra Operations - -```@docs -inner(::MPST, ::MPST) where {MPST <: ITensorMPS.AbstractMPS} -dot(::MPST, ::MPST) where {MPST <: ITensorMPS.AbstractMPS} -loginner(::MPST, ::MPST) where {MPST <: ITensorMPS.AbstractMPS} -logdot(::MPST, ::MPST) where {MPST <: ITensorMPS.AbstractMPS} -inner(::MPS, ::MPO, ::MPS) -dot(::MPS, ::MPO, ::MPS) -inner(::MPO, ::MPS, ::MPO, ::MPS) -dot(::MPO, ::MPS, ::MPO, ::MPS) -norm(::ITensorMPS.AbstractMPS) -normalize(::ITensorMPS.AbstractMPS) -normalize!(::ITensorMPS.AbstractMPS) -lognorm(::ITensorMPS.AbstractMPS) -+(::ITensorMPS.AbstractMPS...) -contract(::MPO, ::MPS) -apply(::MPO, ::MPS) -contract(::MPO, ::MPO) -apply(::MPO, ::MPO) -error_contract(y::MPS, A::MPO, x::MPS) -outer(::MPS, ::MPS) -projector(::MPS) -``` - diff --git a/docs/src/Multithreading.md b/docs/src/Multithreading.md index 6e54bcf86a..4b1f7eb5c6 100644 --- a/docs/src/Multithreading.md +++ b/docs/src/Multithreading.md @@ -88,7 +88,7 @@ Here is a simple example of using block sparse multithreading to speed up a spar tensor contraction: ```julia using BenchmarkTools -using ITensors, ITensorMPS +using ITensors using LinearAlgebra using Strided @@ -103,7 +103,7 @@ function main(; d = 20, order = 4) println() i(n) = Index(QN(0) => d, QN(1) => d; tags = "i$n") - is = IndexSet(i, order ÷ 2) + is = ntuple(i, order ÷ 2) A = random_itensor(is'..., dag(is)...) B = random_itensor(is'..., dag(is)...) diff --git a/docs/src/Observer.md b/docs/src/Observer.md deleted file mode 100644 index 66affcfaa1..0000000000 --- a/docs/src/Observer.md +++ /dev/null @@ -1,200 +0,0 @@ -# [Observer System for DMRG](@id observer) - -An observer is an object which can be passed to the ITensor DMRG -algorithm, to allow measurements to be performed throughout -the DMRG calculation and to set conditions for early stopping -of DMRG. - -The only requirement of an observer is that it is a subtype -of `AbstractObserver`. But to do something interesting, it -should also overload at least one the methods `measure!` -or `checkdone!`. - -A general purpose observer type called [`DMRGObserver`](@ref) is -included with ITensors which already provides some -quite useful features. It accepts a list of strings naming -local operators to be measured at each step of DMRG, with -the results saved for later analysis. It also accepts an -optional energy precision, and stops a DMRG calculation early -if the energy no longer changes to this precision. For more -details about the [`DMRGObserver`](@ref) type, see -the [DMRGObserver](@ref) documentation page. - -## Defining a Custom Observer - -To define a custom observer, just make a struct with -any name and internal fields you would like, and make -this struct a subtype of `AbstractObserver`. - -For example, let's make a type called `DemoObserver` -as: - -```julia -using ITensors, ITensorMPS - -mutable struct DemoObserver <: AbstractObserver - energy_tol::Float64 - last_energy::Float64 - - DemoObserver(energy_tol=0.0) = new(energy_tol,1000.0) -end - -``` - -In this minimal example, our `DemoObserver` -contains a field `energy_tol` which we can use to set -an early-stopping condition for DMRG, and an field -`last_energy` which our observer will use internally -to keep track of changes to the energy after each sweep. - -Now to give our `DemoObserver` type a useful behavior -we need to define overloads of the methods `measure!` -and `checkdone!`. - -### Overloading the `checkdone!` method - -Let's start with the `checkdone!` method. After -each sweep of DMRG, the `checkdone!` method is -passed the observer object, as well as a set of keyword -arguments which currently include: - - energy: the current energy - - psi: the current wavefunction MPS - - sweep: the number of the sweep that just finished - - outputlevel: an integer stating the desired level of output - -If the `checkdone!` function returns `true`, then the DMRG -routine stops (recall that `checkdone!` is called only at the -end of a sweep). - -In our example, we will just compare the `energy` keyword -argument to the `last_energy` variable held inside the `DemoObserver`: - -```julia -function ITensorMPS.checkdone!(o::DemoObserver;kwargs...) - sw = kwargs[:sweep] - energy = kwargs[:energy] - if abs(energy-o.last_energy)/abs(energy) < o.energy_tol - println("Stopping DMRG after sweep $sw") - return true - end - # Otherwise, update last_energy and keep going - o.last_energy = energy - return false -end -``` - -(Recall that in order to properly overload the default behavior, -the `checkdone!` method has to be imported from the ITensors module -or preceded with `ITensors.`) - - -### Overloading the `measure!` method - -The other method that an observer can overload is `measure!`. -This method is called at every step of DMRG, so at every -site and for every sweep. The `measure!` method is passed -the current observer object and a set of keyword arguments -which include: - - energy: the energy after the current step of DMRG - - psi: the current wavefunction MPS - - bond: the bond `b` that was just optimized, corresponding to sites `(b,b+1)` in the two-site DMRG algorithm - - sweep: the current sweep number - - sweep\_is\_done: true if at the end of the current sweep, otherwise false - - half_sweep: the half-sweep number, equal to 1 for a left-to-right, first half sweep, or 2 for the second, right-to-left half sweep - - spec: the Spectrum object returned from factorizing the local superblock wavefunction tensor in two-site DMRG - - outputlevel: an integer specifying the amount of output to show - - projected_operator: projection of the linear operator into the current MPS basis - -For our minimal `DemoObserver` example here, we will just make a `measure!` function -that prints out some of the information above, but in a more realistic setting one -could use the MPS `psi` to perform essentially arbitrary measurements. - -```julia -function ITensorMPS.measure!(o::DemoObserver; kwargs...) - energy = kwargs[:energy] - sweep = kwargs[:sweep] - bond = kwargs[:bond] - outputlevel = kwargs[:outputlevel] - - if outputlevel > 0 - println("Sweep $sweep at bond $bond, the energy is $energy") - end -end -``` - -## Calling DMRG with the Custom Observer - -After defining an observer type and overloading at least one of the -methods `checkdone!` or `measure!` for it, one can construct an -object of this type and pass it to the ITensor [`dmrg`](@ref) function -using the `observer` keyword argument. - -Continuing with our `DemoObserver` example above: - -```julia -obs = DemoObserver(1E-4) # use an energy tolerance of 1E-4 -energy, psi = dmrg(H,psi0,sweeps; observer=obs, outputlevel=1) -``` - -## Complete Sample Code - -```julia -using ITensors, ITensorMPS - -mutable struct DemoObserver <: AbstractObserver - energy_tol::Float64 - last_energy::Float64 - - DemoObserver(energy_tol=0.0) = new(energy_tol,1000.0) -end - -function ITensorMPS.checkdone!(o::DemoObserver;kwargs...) - sw = kwargs[:sweep] - energy = kwargs[:energy] - if abs(energy-o.last_energy)/abs(energy) < o.energy_tol - println("Stopping DMRG after sweep $sw") - return true - end - # Otherwise, update last_energy and keep going - o.last_energy = energy - return false -end - -function ITensorMPS.measure!(o::DemoObserver; kwargs...) - energy = kwargs[:energy] - sweep = kwargs[:sweep] - bond = kwargs[:bond] - outputlevel = kwargs[:outputlevel] - - if outputlevel > 0 - println("Sweep $sweep at bond $bond, the energy is $energy") - end -end - -let - N = 10 - etol = 1E-4 - - s = siteinds("S=1/2",N) - - a = OpSum() - for n=1:N-1 - a += "Sz",n,"Sz",n+1 - a += 0.5,"S+",n,"S-",n+1 - a += 0.5,"S-",n,"S+",n+1 - end - H = MPO(a,s) - psi0 = random_mps(s;linkdims=4) - - nsweeps = 5 - cutoff = 1E-8 - maxdim = [10,20,100] - - obs = DemoObserver(etol) - - println("Starting DMRG") - energy, psi = dmrg(H,psi0; nsweeps, cutoff, maxdim, observer=obs, outputlevel=1) - - return -end -``` diff --git a/docs/src/OpSum.md b/docs/src/OpSum.md deleted file mode 100644 index 89234e2b54..0000000000 --- a/docs/src/OpSum.md +++ /dev/null @@ -1,14 +0,0 @@ -# OpSum - -## Description - -```@docs -OpSum -``` - -## Methods - -```@docs -add! -MPO(::OpSum,::Vector{<:Index}) -``` diff --git a/docs/src/ProjMPO.md b/docs/src/ProjMPO.md deleted file mode 100644 index 8010775b3f..0000000000 --- a/docs/src/ProjMPO.md +++ /dev/null @@ -1,23 +0,0 @@ -# ProjMPO - -## Description - -```@docs -ProjMPO -``` - -## Methods - -```@docs -product(::ProjMPO,::ITensor) -position!(::ProjMPO, ::MPS, ::Int) -noiseterm(::ProjMPO,::ITensor,::String) -``` - -## Properties - -```@docs -length(::ProjMPO) -eltype(::ProjMPO) -size(::ProjMPO) -``` diff --git a/docs/src/ProjMPOSum.md b/docs/src/ProjMPOSum.md deleted file mode 100644 index 8a713b269e..0000000000 --- a/docs/src/ProjMPOSum.md +++ /dev/null @@ -1,22 +0,0 @@ -# ProjMPOSum - -## Description - -```@docs -ProjMPOSum -``` - -## Methods - -```@docs -product(::ProjMPOSum,::ITensor) -position!(::ProjMPOSum, ::MPS, ::Int) -noiseterm(::ProjMPOSum,::ITensor,::String) -``` - -## Properties - -```@docs -eltype(::ProjMPOSum) -size(::ProjMPOSum) -``` diff --git a/docs/src/QNTricks.md b/docs/src/QNTricks.md deleted file mode 100644 index abf7a8239e..0000000000 --- a/docs/src/QNTricks.md +++ /dev/null @@ -1,81 +0,0 @@ -# Symmetric (QN Conserving) Tensors: Background and Usage - -Here is a collection of background material and example codes for understanding how symmetric tensors (tensors with conserved quantum numbers) work in ITensors.jl - -## Combiners and Symmetric Tensors - -In ITensors.jl, combiners are special sparse tensors that represent the action of taking the tensor product of one or more indices. It generalizes the idea of reshaping and permuting. For dense ITensors, a combiner is just the action of permuting and reshaping the data of the tensor. For symmetric tensors (quantum number conserving tensors represented as block sparse tensors), the combiner also fuses symmetry sectors together. They can be used for various purposes. Generally they are used internally in the library, for example in order to reshape a high order ITensor into an order 2 ITensor to perform a matrix decomposition like an SVD or eigendecomposition. - -For example: -```@repl -using ITensors - -# This is a short code showing how a combiner -# can be used to "flip" the direction of an Index -i = Index([QN(0) => 2, QN(1) => 3], "i") -j = Index([QN(0) => 2, QN(1) => 3], "j") -A = random_itensor(i, dag(j)) -C = combiner(i, dag(j); tags = "c", dir = dir(i)) -inds(A) -inds(A * C) -``` -You can see that the combiner reshapes the indices of `A` into a single Index that contains the tensor product of the two input spaces. The spaces have size `QN(-1) => 2 * 3`, `QN(0) => 2 * 2 + 3 * 3`, and `QN(0) => 2 * 3` (determined from all of the combinations of combining the sectors of the different indices, where the QNs are added and the block dimensions are multiplied). The ordering of the sectors is determined internally by ITensors.jl. - -You can also use a combiner on a single Index, which can be helpful for changing the direction of an Index or combining multiple sectors of the same symmetry into a single sector: -```@repl -using ITensors - -# This is a short code showing how a combiner -# can be used to "flip" the direction of an Index -i = Index([QN(0) => 2, QN(1) => 3], "i") -j = dag(Index([QN(0) => 2, QN(1) => 3], "j")) -A = random_itensor(i, j) -C = combiner(j; tags = "jflip", dir = -dir(j)) -inds(A) -inds(A * C) -``` -Unless you are writing very specialized custom code with symmetric tensors, this is generally not needed. - -## Block Sparsity and Quantum Numbers - -In general, not all blocks that are allowed according to the flux will actually exist in the tensor (which helps in many cases for efficiency). Usually this would happen when the tensor is first constructed and not all blocks are explicitly set: -```@repl -using ITensors - -i = Index([QN(0) => 1, QN(1) => 1]) -A = ITensor(i', dag(i)); -A[2, 2] = 1.0; -@show A; -D, U = eigen(A; ishermitian=true); -@show D; -@show U; -``` -If we had set `A[1, 1] = 0.0` as well, then all of the allowed blocks (according to the flux `QN(0)` would exist and would be included in the eigendecomposition: -```@repl -using ITensors - -i = Index([QN(0) => 1, QN(1) => 1]) -A = ITensor(i', dag(i)); -A[2, 2] = 1.0; -A[1, 1] = 0.0; -@show A; -D, U = eigen(A; ishermitian=true); -@show D; -@show U; -``` -"Missing" blocks can also occur with tensor contractions, since the final blocks of the output tensor are made from combinations of contractions of blocks from the input tensors, and there is no guarantee that all flux-consistent blocks will end up in the result: -```@repl -using ITensors - -i = Index([QN(0) => 1, QN(1) => 1]) -j = Index([QN(0) => 1]) -A = ITensor(i, dag(j)); -A[2, 1] = 1.0; -@show A; -A2 = prime(A, i) * dag(A); -@show A2; -D, U = eigen(A2; ishermitian=true); -@show D; -@show U; -``` - diff --git a/docs/src/SiteType.md b/docs/src/SiteType.md deleted file mode 100644 index 42f4d116e1..0000000000 --- a/docs/src/SiteType.md +++ /dev/null @@ -1,17 +0,0 @@ -# SiteType and op, state, val functions - -## Description - -```@docs -SiteType -``` - -## Methods - -```@docs -op -state -val -space -``` - diff --git a/docs/src/Sweeps.md b/docs/src/Sweeps.md deleted file mode 100644 index 889caa7a2e..0000000000 --- a/docs/src/Sweeps.md +++ /dev/null @@ -1,26 +0,0 @@ -# Sweeps - -```@docs -Sweeps -Sweeps(nsw::Int, d::AbstractMatrix) -``` - -## Modifying Sweeps Objects - -```@docs -setmaxdim! -setcutoff! -setnoise! -setmindim! -``` - -## Getting Sweeps Object Data - -```@docs -nsweep(sw::Sweeps) -maxdim(sw::Sweeps,n::Int) -cutoff(sw::Sweeps,n::Int) -noise(sw::Sweeps,n::Int) -mindim(sw::Sweeps,n::Int) -``` - diff --git a/docs/src/examples/DMRG.md b/docs/src/examples/DMRG.md deleted file mode 100644 index 48d51feb11..0000000000 --- a/docs/src/examples/DMRG.md +++ /dev/null @@ -1,552 +0,0 @@ -# DMRG Code Examples - -## Perform a basic DMRG calculation - -Because tensor indices in ITensor have unique identities, before we can make a Hamiltonian -or a wavefunction we need to construct a "site set" which will hold the site indices defining -the physical Hilbert space: - -```julia -using ITensors, ITensorMPS -N = 100 -sites = siteinds("S=1",N) -``` - -Here we have chosen to create a Hilbert space of N spin 1 sites. The string "S=1" -denotes a special Index tag which hooks into a system that knows "S=1" indices have -a dimension of 3 and how to create common physics operators like "Sz" for them. - -Next we'll make our Hamiltonian matrix product operator (MPO). A very -convenient way to do this is to use the OpSum helper type which lets -us input a Hamiltonian (or any sum of local operators) in similar notation -to pencil-and-paper notation: - -```julia -os = OpSum() -for j=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,sites) -``` - -In the last line above we convert the OpSum helper object to an actual MPO. - -Before beginning the calculation, we need to specify how many DMRG sweeps to do and -what schedule we would like for the parameters controlling the accuracy. -These parameters can be specified as follows: - -```julia -nsweeps = 5 # number of sweeps is 5 -maxdim = [10,20,100,100,200] # gradually increase states kept -cutoff = [1E-10] # desired truncation error -``` - -The random starting wavefunction `psi0` must be defined in the same Hilbert space -as the Hamiltonian, so we construct it using the same collection of site indices: - -```julia -psi0 = random_mps(sites;linkdims=2) -``` - -Here we have made a random MPS of bond dimension 2. We could have used a random product -state instead, but choosing a slightly larger bond dimension can help DMRG avoid getting -stuck in local minima. We could also set psi to some specific initial state using the -`MPS` constructor, which is actually required if we were conserving QNs. - -Finally, we are ready to call DMRG: - -```julia -energy,psi = dmrg(H,psi0;nsweeps,maxdim,cutoff) -``` - -When the algorithm is done, it returns the ground state energy as the variable `energy` and an MPS -approximation to the ground state as the variable `psi`. - -Below you can find a complete working code that includes all of these steps: - -```julia -using ITensors, ITensorMPS - -let - N = 100 - sites = siteinds("S=1",N) - - os = OpSum() - for j=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,sites) - - nsweeps = 5 # number of sweeps is 5 - maxdim = [10,20,100,100,200] # gradually increase states kept - cutoff = [1E-10] # desired truncation error - - psi0 = random_mps(sites;linkdims=2) - - energy,psi = dmrg(H,psi0;nsweeps,maxdim,cutoff) - - return -end -``` - -## Using a Custom Observer for DMRG - -An Observer is any object which can be used to perform custom measurements throughout -a DMRG calculation and to stop a DMRG calculation early. Because an Observer has -access to the entire wavefunction at every step, a wide range of customization is -possible. - -For detailed examples of making custom Observers, see the [Observer](@ref observer) -section of the documentation. - - -## DMRG Calculation with Mixed Local Hilbert Space Types - -The following fully-working example shows how to set up a calculation -mixing S=1/2 and S=1 spins on every other site of a 1D system. The -Hamiltonian involves Heisenberg spin interactions with adjustable -couplings between sites of the same spin or different spin. - -Note that the only difference from a regular ITensor DMRG calculation -is that the `sites` array has Index objects which alternate in dimension -and in which physical tag type they carry, whether `"S=1/2"` or `"S=1"`. -(Try printing out the sites array to see!) -These tags tell the OpSum system which local operators to use for these -sites when building the Hamiltonian MPO. - -```julia -using ITensors, ITensorMPS - -let - N = 100 - - # Make an array of N Index objects with alternating - # "S=1/2" and "S=1" tags on odd versus even sites - # (The first argument n->isodd(n) ... is an - # on-the-fly function mapping integers to strings) - sites = siteinds(n->isodd(n) ? "S=1/2" : "S=1",N) - - # Couplings between spin-half and - # spin-one sites: - Jho = 1.0 # half-one coupling - Jhh = 0.5 # half-half coupling - Joo = 0.5 # one-one coupling - - os = OpSum() - for j=1:N-1 - os += 0.5*Jho,"S+",j,"S-",j+1 - os += 0.5*Jho,"S-",j,"S+",j+1 - os += Jho,"Sz",j,"Sz",j+1 - end - for j=1:2:N-2 - os += 0.5*Jhh,"S+",j,"S-",j+2 - os += 0.5*Jhh,"S-",j,"S+",j+2 - os += Jhh,"Sz",j,"Sz",j+2 - end - for j=2:2:N-2 - os += 0.5*Joo,"S+",j,"S-",j+2 - os += 0.5*Joo,"S-",j,"S+",j+2 - os += Joo,"Sz",j,"Sz",j+2 - end - H = MPO(os,sites) - - nsweeps = 10 - maxdim = [10,10,20,40,80,100,140,180,200] - cutoff = [1E-8] - - psi0 = random_mps(sites;linkdims=4) - - energy,psi = dmrg(H,psi0;nsweeps,maxdim,cutoff) - - return -end -``` - -## Use a Sum of MPOs in DMRG - -One version of the ITensor `dmrg` function accepts an array of MPOs -`[H1,H2,H3]` (or any number of MPOs you want). This version of DMRG -will find the ground state of `H1+H2+H3`. Internally it does not -actually sum these MPOs, but loops over them during each step of -the "eigensolver" at the core of the DMRG algorithm, so it is usually -more efficient than if the MPOs had been summed together into a single MPO. - -To use this version of DMRG, say you have MPOs `H1`, `H2`, and `H3`. -Then call DMRG like this: -```julia -energy,psi = dmrg([H1,H2,H3],psi0;nsweeps,maxdim,cutoff) -``` - -## Make a 2D Hamiltonian for DMRG - -You can use the OpSum system to make 2D Hamiltonians -much in the same way you make 1D Hamiltonians: by looping over -all of the bonds and adding the interactions on these bonds to -the OpSum. - -To help with the logic of 2D lattices, ITensor pre-defines -some helper functions which -return an array of bonds. Each bond object has an -"s1" field and an "s2" field which are the integers numbering -the two sites the bond connects. -(You can view the source for these functions at [this link](https://github.com/ITensor/ITensors.jl/blob/main/src/lib/ITensorMPS/src/lattices/lattices.jl).) - -The two provided functions currently are `square_lattice` and -`triangular_lattice`. It is not hard to write your own similar lattice -functions as all they have to do is define an array of `ITensors.LatticeBond` -structs or even a custom struct type you wish to define. We welcome any -user contributions of other lattices that ITensor does not currently offer. - -Each lattice function takes an optional named argument -"yperiodic" which lets you request that the lattice should -have periodic boundary conditions around the y direction, making -the geometry a cylinder. - -**Full example code:** - -```julia -using ITensors, ITensorMPS - -let - Ny = 6 - Nx = 12 - - N = Nx*Ny - - sites = siteinds("S=1/2", N; - conserve_qns = true) - - # Obtain an array of LatticeBond structs - # which define nearest-neighbor site pairs - # on the 2D square lattice (wrapped on a cylinder) - lattice = square_lattice(Nx, Ny; yperiodic = false) - - # Define the Heisenberg spin Hamiltonian on this lattice - os = OpSum() - for b in lattice - os += 0.5, "S+", b.s1, "S-", b.s2 - os += 0.5, "S-", b.s1, "S+", b.s2 - os += "Sz", b.s1, "Sz", b.s2 - end - H = MPO(os,sites) - - state = [isodd(n) ? "Up" : "Dn" for n=1:N] - # Initialize wavefunction to a random MPS - # of bond-dimension 10 with same quantum - # numbers as `state` - psi0 = random_mps(sites,state;linkdims=20) - - nsweeps = 10 - maxdim = [20,60,100,100,200,400,800] - cutoff = [1E-8] - - energy,psi = dmrg(H,psi0;nsweeps,maxdim,cutoff) - - return -end -``` - -## Compute excited states with DMRG - -ITensor DMRG accepts additional MPS wavefunctions as a optional, extra argument. -These additional 'penalty states' are provided as an array of MPS just -after the Hamiltonian, like this: - -```julia -energy,psi3 = dmrg(H,[psi0,psi1,psi2],psi3_init;nsweeps,maxdim,cutoff) -``` - -Here the penalty states are `[psi0,psi1,psi2]`. -When these are provided, the DMRG code minimizes the -energy of the current MPS while also reducing its overlap -(inner product) with the previously provided MPS. If these overlaps become sufficiently small, -then the computed MPS is an excited state. So by finding the ground -state, then providing it to DMRG as a "penalty state" or previous state -one can compute the first excited state. Then providing both of these, one can -get the second excited state, etc. - -A keyword argument called `weight` can also be provided to -the `dmrg` function when penalizing overlaps to previous states. The -`weight` parameter is multiplied by the overlap with the previous states, -so sets the size of the penalty. It should be chosen at least as large -as the (estimated) gap between the ground and first excited states. -Otherwise the optimal value of the weight parameter is not so obvious, -and it is best to try various weights during initial test calculations. - -Note that when the system has conserved quantum numbers, a superior way -to find excited states can be to find ground states of quantum number (or symmetry) -sectors other than the one containing the absolute ground state. In that -context, the penalty method used below is a way to find higher excited states -within the same quantum number sector. - -**Full Example code:** - -```julia -using ITensors, ITensorMPS - -let - N = 20 - - sites = siteinds("S=1/2",N) - - h = 4.0 - - weight = 20*h # use a large weight - # since gap is expected to be large - - - # - # Use the OpSum feature to create the - # transverse field Ising model - # - # Factors of 4 and 2 are to rescale - # spin operators into Pauli matrices - # - os = OpSum() - for j=1:N-1 - os -= 4,"Sz",j,"Sz",j+1 - end - for j=1:N - os -= 2*h,"Sx",j; - end - H = MPO(os,sites) - - - # - # Make sure to do lots of sweeps - # when finding excited states - # - nsweeps = 30 - maxdim = [10,10,10,20,20,40,80,100,200,200] - cutoff = [1E-8] - noise = [1E-6] - - # - # Compute the ground state psi0 - # - psi0_init = random_mps(sites;linkdims=2) - energy0,psi0 = dmrg(H,psi0_init;nsweeps,maxdim,cutoff,noise) - - println() - - # - # Compute the first excited state psi1 - # - psi1_init = random_mps(sites;linkdims=2) - energy1,psi1 = dmrg(H,[psi0],psi1_init;nsweeps,maxdim,cutoff,noise,weight) - - # Check psi1 is orthogonal to psi0 - @show inner(psi1,psi0) - - - # - # The expected gap of the transverse field Ising - # model is given by Eg = 2*|h-1| - # - # (The DMRG gap will have finite-size corrections) - # - println("DMRG energy gap = ",energy1-energy0); - println("Theoretical gap = ",2*abs(h-1)); - - println() - - # - # Compute the second excited state psi2 - # - psi2_init = random_mps(sites;linkdims=2) - energy2,psi2 = dmrg(H,[psi0,psi1],psi2_init;nsweeps,maxdim,cutoff,noise,weight) - - # Check psi2 is orthogonal to psi0 and psi1 - @show inner(psi2,psi0) - @show inner(psi2,psi1) - - return -end -``` - -## Printing the Entanglement Entropy at Each Step - -To obtain the entanglement entropy of an MPS at each step during a DMRG calculation, -you can use the [Observer](@ref observer) system to make a custom observer object that prints out -this information. - -First we define our custom observer type, `EntanglementObserver`, and overload the `measure!` function -for it: - -```julia -using ITensors, ITensorMPS - -mutable struct EntanglementObserver <: AbstractObserver -end - -function ITensorMPS.measure!(o::EntanglementObserver; bond, psi, half_sweep, kwargs...) - wf_center, other = half_sweep==1 ? (psi[bond+1],psi[bond]) : (psi[bond],psi[bond+1]) - U,S,V = svd(wf_center, uniqueinds(wf_center,other)) - SvN = 0.0 - for n=1:dim(S, 1) - p = S[n,n]^2 - SvN -= p * log(p) - end - println(" Entanglement across bond $bond = $SvN") -end -``` - -The `measure!` function grabs certain helpful keywords passed to it by DMRG, such as what bond DMRG -has just finished optimizing. - -Here is a complete sample code including constructing the observer and passing it to DMRG: - -```julia -using ITensors, ITensorMPS - -mutable struct EntanglementObserver <: AbstractObserver -end - -function ITensorMPS.measure!(o::EntanglementObserver; bond, psi, half_sweep, kwargs...) - wf_center, other = half_sweep==1 ? (psi[bond+1],psi[bond]) : (psi[bond],psi[bond+1]) - U,S,V = svd(wf_center, uniqueinds(wf_center,other)) - SvN = 0.0 - for n=1:dim(S, 1) - p = S[n,n]^2 - SvN -= p * log(p) - end - println(" Entanglement across bond $bond = $SvN") -end - -let - N = 100 - - s = siteinds("S=1/2",N) - - a = OpSum() - for n=1:N-1 - a += "Sz",n,"Sz",n+1 - a += 0.5,"S+",n,"S-",n+1 - a += 0.5,"S-",n,"S+",n+1 - end - H = MPO(a,s) - psi0 = random_mps(s;linkdims=4) - - nsweeps = 5 - maxdim = [10,20,80,160] - cutoff = 1E-8 - - observer = EntanglementObserver() - - energy,psi = dmrg(H,psi0;nsweeps,maxdim,cutoff,observer,outputlevel=2) - - return -end -``` - -Example output: -``` -... -Sweep 2, half 2, bond (35,36) energy=-44.08644657103751 - Truncated using cutoff=1.0E-08 maxdim=20 mindim=1 - Trunc. err=2.54E-07, bond dimension 20 - Entanglement across bond 35 = 0.7775882479059774 -Sweep 2, half 2, bond (34,35) energy=-44.086696891668424 - Truncated using cutoff=1.0E-08 maxdim=20 mindim=1 - Trunc. err=2.12E-07, bond dimension 20 - Entanglement across bond 34 = 0.7103532704635472 -Sweep 2, half 2, bond (33,34) energy=-44.08696190368391 - Truncated using cutoff=1.0E-08 maxdim=20 mindim=1 - Trunc. err=1.29E-07, bond dimension 20 - Entanglement across bond 33 = 0.7798362911744212 -... -``` - -If you only want to see the maximum entanglement during each sweep, you can add a field to the EntanglementObserver -object that saves the maximum value encountered so far and keep overwriting this field, printing out the most recently observed maximum at the end of each sweep. - - -## Monitoring the Memory Usage of DMRG - -To monitor how much memory (RAM) a DMRG calculation is using while it is running, -you can use the [Observer](@ref observer) system to make a custom observer object that prints out -this information. Also the `Base.summarysize` function, which returns the size -in bytes of any Julia object is very helpful here. - -First we define our custom observer type, `SizeObserver`, and overload the `measure!` function -for it: - -```julia -using ITensors, ITensorMPS - -mutable struct SizeObserver <: AbstractObserver -end - -function ITensorMPS.measure!(o::SizeObserver; bond, half_sweep, psi, projected_operator, kwargs...) - if bond==1 && half_sweep==2 - psi_size = Base.format_bytes(Base.summarysize(psi)) - PH_size = Base.format_bytes(Base.summarysize(projected_operator)) - println("|psi| = $psi_size, |PH| = $PH_size") - end -end -``` - -The `measure!` function grabs certain helpful keywords passed to it by DMRG, checking -`if bond==1 && half_sweep==2` so that it only runs when at the end of a full sweep. - -When it runs, it calls `Base.summarysize` on the wavefunction `psi` object and the `projected_operator` object. The `projected_operator`, which is the matrix (Hamiltonian) wrapped into the current MPS basis, is usually the largest-sized object in a DMRG calculation. The code also uses `Base.format_bytes` to turn an integer representing bytes into a human-readable string. - -Here is a complete sample code including constructing the observer and passing it to DMRG: - -```julia -using ITensors, ITensorMPS - -mutable struct SizeObserver <: AbstractObserver -end - -function ITensorMPS.measure!(o::SizeObserver; bond, sweep, half_sweep, psi, projected_operator, kwargs...) - if bond==1 && half_sweep==2 - psi_size = Base.format_bytes(Base.summarysize(psi)) - PH_size = Base.format_bytes(Base.summarysize(projected_operator)) - println("After sweep $sweep, |psi| = $psi_size, |PH| = $PH_size") - end -end - -let - N = 100 - - s = siteinds("S=1/2",N) - - a = OpSum() - for n=1:N-1 - a += "Sz",n,"Sz",n+1 - a += 0.5,"S+",n,"S-",n+1 - a += 0.5,"S-",n,"S+",n+1 - end - H = MPO(a,s) - psi0 = random_mps(s;linkdims=4) - - nsweeps = 5 - maxdim = [10,20,80,160] - cutoff = 1E-8 - - obs = SizeObserver() - - energy,psi = dmrg(H,psi0;nsweeps,maxdim,cutoff,observer=obs) - - return -end -``` - -Example output: -``` -After sweep 1, |psi| = 211.312 KiB, |PH| = 593.984 KiB -After sweep 1 energy=-43.95323393592883 maxlinkdim=10 maxerr=8.26E-06 time=0.098 -After sweep 2, |psi| = 641.000 KiB, |PH| = 1.632 MiB -After sweep 2 energy=-44.10791340895817 maxlinkdim=20 maxerr=7.39E-07 time=0.132 -After sweep 3, |psi| = 1.980 MiB, |PH| = 5.066 MiB -After sweep 3 energy=-44.12593605906466 maxlinkdim=44 maxerr=9.96E-09 time=0.256 -After sweep 4, |psi| = 2.863 MiB, |PH| = 7.246 MiB -After sweep 4 energy=-44.127710946536645 maxlinkdim=56 maxerr=9.99E-09 time=0.445 -After sweep 5, |psi| = 3.108 MiB, |PH| = 7.845 MiB -After sweep 5 energy=-44.127736798226536 maxlinkdim=57 maxerr=9.98E-09 time=0.564 -``` diff --git a/docs/src/examples/MPSandMPO.md b/docs/src/examples/MPSandMPO.md deleted file mode 100644 index c8ee559261..0000000000 --- a/docs/src/examples/MPSandMPO.md +++ /dev/null @@ -1,454 +0,0 @@ -# MPS and MPO Examples - -The following examples demonstrate operations available in ITensor -to work with [matrix product state (MPS)](http://tensornetwork.org/mps/) -(or tensor train) and matrix product operator (MPO) tensor networks. - -## Creating an MPS from a Tensor - -![](mps_from_tensor.png) - -A matrix product state (MPS) made of N tensors, each with -one site or physical index, is a way of representing a single -tensor with N indices. One way of obtaining the MPS form of an -N-index tensor `T` is by repeatedly factorizing `T` into N -separate tensors using a factorization such as the [Singular Value Decomposition](@ref) (SVD). -This algorithm for obtaining an MPS is known in the mathematics -literature as the "tensor train SVD" or "TT-SVD" algorithm. - -To turn an N-index (order-N) tensor T into an MPS, you can just -construct an MPS by passing T as the first argument, along with -keyword arguments that control the approximations used in factorizing -T. Let's look at a few specific cases. - -#### ITensor to MPS Example - -If you have a tensor `T` which is an ITensor and has indices `i,j,k,l,m`, -you can create an MPS approximation of `T` where the MPS has site indices -`i,j,k,l,m` as follows: - -```julia -cutoff = 1E-8 -maxdim = 10 -T = random_itensor(i,j,k,l,m) -M = MPS(T,(i,j,k,l,m);cutoff=cutoff,maxdim=maxdim) -``` - -Here we used a random ITensor for illustrative purposes, but it could be any ITensor and -typically tensors with additional structure are more well approximated by MPS. - -#### Julia Tensor to MPS Example - -Another situation could be where you have a Julia array or Julia tensor of -dimension ``d^N`` and want to approximate it as an MPS with ``N`` site indices, -each of dimension ``d``. For example, we could have the following random Julia -array of dimension ``2\times 2\times 2 \times 2 \times 2``: - -```julia -d = 2 -N = 5 -A = randn(d,d,d,d,d) -``` - -Alternatively, the array could be just a one dimensional array of length ``d^N``: - -```julia -A = randn(d^N) -``` - -To convert this array to an MPS, we will first need a collection of Index objects -to use as the site indices of the MPS. We can conveniently construct an array of -four indices of dimension 2 as follows: - -```julia -sites = siteinds(d,N) -``` - -Finally, we can pass our array `A` and our `sites` to the MPS constructor along with -parameters controlling the truncation level of the factorizations used: - -```julia -cutoff = 1E-8 -maxdim = 10 -M = MPS(A,sites;cutoff=cutoff,maxdim=maxdim) -``` - -## Obtaining Elements of a Tensor Represented by an MPS - -A matrix product state (MPS) or tensor train (TT) is a format for representing a large tensor having N indices in terms of N smaller tensors. Given an MPS represeting a tensor T -we can obtain a particular element ``T^{s_1 s_2 s_3 \cdots s_N}`` -of that tensor using code similar to the following code below. - -In the example code below we will obtain the element ``T^{1,2,1,1,2,1,2,2,2,1}`` of the tensor T -which is (implicitly) defined by the MPS psi: - -```@example mps_element -using ITensors, ITensorMPS -let # hide -N = 10 -s = siteinds(2,N) -chi = 4 -psi = random_mps(s;linkdims=chi) - -# Make an array of integers of the element we -# want to obtain -el = [1,2,1,1,2,1,2,2,2,1] - -V = ITensor(1.) -for j=1:N - V *= (psi[j]*state(s[j],el[j])) -end -v = scalar(V) - -# v is the element we wanted to obtain: -@show v -end # hide -``` - -The call to `state(s[j],el[j])` in the code above makes a single-index ITensor -with the Index `s[j]` and the entry at location `el[j]` set to 1.0, with all other -entries set to 0.0. Contracting this tensor with the MPS tensor at site `j` -can be viewed as "clamping" or "fixing" the index to a set value. The resulting -tensors are contracted sequentially, overwriting the ITensor `V`, and the final -scalar value of `V` is the tensor element we seek. - -See below for a visual depiction of what the above code is doing: - -![](mps_element.png) - -## Expected Value of Local Operators - -When using an MPS to represent a quantum wavefunction ``|\psi\rangle`` -a common operation is computing the expected value ``\langle\psi|\hat{A}_j|\psi\rangle`` -of a local operator ``\hat{A}_j`` acting on site ``j``. This can be accomplished -efficiently and conveniently using the [`expect`](@ref) function as: - -```julia -Avals = expect(psi,"A") -``` - -where `"A"` must be an operator associated with the physical site type, or site tags, of -the sites of the MPS `psi`. For example, the operator name could be -`"Sz"` for spin sites or `"Ntot"` for electron sites. -(For more information about defining such operators yourself, -see the section on [Extending op Function Definitions](@ref custom_op).) - -As a concrete example, consider computing the expectation value of ``S^z_j`` on -every site of an MPS representing a system of N spins of size ``S=1/2``. In the -following example we will use a random MPS of bond dimension ``\chi=4`` but the -MPS could be obtained other ways such as through a DMRG calculation. - -```@example expect -using ITensors, ITensorMPS -N = 10 -chi = 4 -sites = siteinds("S=1/2",N) -psi = random_mps(sites;linkdims=chi) -magz = expect(psi,"Sz") -for (j,mz) in enumerate(magz) - println("$j $mz") -end -``` - -![](mps_expect.png) - -## Expected Values of MPO Operators - -When using an MPS to represent a quantum wavefunction ``|\psi\rangle`` -another common operation is computing the expected value ``\langle\psi|W|\psi\rangle`` -of an operator ``W`` which is represented as a matrix product operator (MPO) tensor network. -A key example could be the Hamiltonian defining a quantum system. - -Given an MPO `W` and an MPS `psi`, you can compute ``\langle\psi|W|\psi\rangle`` -by using the function `inner` as follows: -```julia -ex_W = inner(psi',W,psi) -``` -which will return a scalar that may be either real or complex, depending on the properties of -`psi` and `W`. - -## Computing Correlation Functions - -In addition to expected values of local operators -discussed above, another type of observable that is very important -in physics studies are correlation functions of the form - -```math -C_{ij} = \langle\psi| A_i B_j |\psi\rangle -``` - -These can be computed efficiently for an MPS `psi` in ITensor -using the [`correlation_matrix`](@ref) function: - -```julia -C = correlation_matrix(psi,"A","B") -``` - -where `"A"` and `"B"` must be an operator names associated with the physical site type, -or site tags, of the sites of the MPS `psi`. For example, these strings could be -`"Sz"`, `"S+"`, or `"S-"` for spin sites, or `"Cdagup"` and `"Cup"` for electron sites. -(For more information about defining such operators yourself, -see the section on [Extending op Function Definitions](@ref custom_op).) - -As a concrete example, say we have an MPS `psi` for a system of spins and -want to compute the correlator ``\langle\psi|S^z_i S^z_j|\psi\rangle``. -We can compute this as: - -```julia -zzcorr = correlation_matrix(psi,"Sz","Sz") -``` - -![](mps_zz_correlation.png) - -See the [`correlation_matrix`](@ref) docs for more details about additional arguments you can pass -to this function. - -## Applying a Single-site Operator to an MPS - -In many applications one needs to modify a matrix product -state (MPS) by multiplying it with an operator that acts -only on a single site. This is actually a very straightforward -operation and this formula shows you how to do it in ITensor. - -Say we have an operator ``G^{s'_3}_{s_3}`` which acts non-trivially on site 3 of our MPS `psi` -as in the following diagram: - -![](mps_onesite_figures/operator_app_mps.png) - -To carry out this operation, contract the operator G with the MPS tensor for site 3, -removing the prime from the ``s'_3`` index afterward: - -![](mps_onesite_figures/operator_contract.png) - -```julia -newA = G * psi[3] -newA = noprime(newA) -``` - -Finally, put the new tensor back into MPS `psi` to update its third MPS tensor: - -```julia -psi[3] = newA -``` - -Afterward, we can visualize the modified MPS as: - -![](mps_onesite_figures/updated_mps.png) - -As a technical note, if you are working in a context where gauge or orthogonality -properties of the MPS are important, such as in time evolution using two-site gates, -then you may want to call `psi = orthogonalize(psi, 3)` -before modifying the tensor at site 3, which will ensure that the MPS remains in a -well-defined orthogonal gauge centered on site 3. Modifying a tensor which is left- or right-orthogonal -(i.e. not the "center" tensor of the gauge) will destroy the gauge condition and -require extra operations to restore it. (Calling `orthogonalize` method will automatically -fix this but will have to do extra work to do so.) - - -## Applying a Two-site Operator to an MPS - -A very common operation with matrix product states (MPS) is -multiplication by a two-site operator or "gate" which modifies -the MPS. This procedure can be carried out in an efficient, -controlled way which is adaptive in the MPS bond dimension. - -Say we have an operator ``G^{s'_3 s'_4}_{s_3 s_4}`` which -is our gate and which acts on physical sites 3 and 4 of our MPS `psi`, -as in the following diagram: - -![](twosite_figures/gate_app_mps.png) - -To apply this gate in a controlled manner, first 'gauge' the MPS `psi` such -that either site 3 or 4 is the *orthogonality center*. Here we make site 3 -the center: - -```julia -psi = orthogonalize(psi, 3) -``` - -![](twosite_figures/gate_gauge.png) - -The other MPS tensors are now either left-orthogonal or right-orthogonal and can be -left out of further steps without producing incorrect results. - -Next, contract the gate tensor G with the MPS tensors for sites 3 and 4 - -![](twosite_figures/gate_contract.png) - -```julia -wf = (psi[3] * psi[4]) * G -wf = noprime(wf) -``` - -Finally, use the singular value decomposition (SVD) to factorize the -resulting tensor, multiplying the singular values into either U or V. -Assign these two tensors back into the MPS to update it. - -![](twosite_figures/gate_svd.png) - -```julia -inds3 = uniqueinds(psi[3],psi[4]) -U,S,V = svd(wf,inds3,cutoff=1E-8) -psi[3] = U -psi[4] = S*V -``` - -The call to `uniqueinds(psi[3])` analyzes the indices of `psi[3]` and `psi[4]` -and finds any which are unique to just `psi[3]`, saving this collection of indices as `inds3`. -Passing this collection of indices to the `svd` function tells it to treat any indices -that are unique to `psi[3]` as the indices which should go onto the `U` tensor afterward. -We also set a truncation error cutoff of 1E-8 in the call to `svd` to truncate -the smallest singular values and control the size of the resulting MPS. -Other cutoff values can be used, depending on the desired accuracy, -as well as limits on the maximum bond dimension (`maxdim` keyword argument). - -**Complete code example** - -```julia -using ITensors, ITensorMPS - -psi = orthogonalize(psi, 3) - -wf = (psi[3] * psi[4]) * G -wf = noprime(wf) - -inds3 = uniqueinds(psi[3], psi[4]) -U, S, V = svd(wf, inds3; cutoff=1E-8) -psi[3] = U -psi[4] = S * V -``` - -## Computing the Entanglement Entropy of an MPS - -A key advantage of using the matrix product state (MPS) format to represent quantum wavefunctions is that it allows one to efficiently compute the entanglement entropy of any left-right bipartition of the system in one dimension, or for a two-dimensional system any "cut" along the MPS path. - -Say that we have obtained an MPS `psi` of length N and we wish to compute the entanglement entropy of a bipartition of the system into a region "A" which consists of sites 1,2,...,b and a region B consisting of sites b+1,b+2,...,N. - -Then the following code formula can be used to accomplish this task: - -```julia -psi = orthogonalize(psi, b) -U,S,V = svd(psi[b], (linkinds(psi, b-1)..., siteinds(psi, b)...)) -SvN = 0.0 -for n=1:dim(S, 1) - p = S[n,n]^2 - SvN -= p * log(p) -end -``` - -As a brief explanation of the code above, the call to `psi = orthogonalize(psi, b)` -shifts the orthogonality center to site `b` of the MPS. - -The call to the `svd` routine says to treat the link (virtual or bond) Index connecting the b'th MPS tensor `psi[b]` and the b'th physical Index as "row" indices for the purposes of the SVD (these indices will end up on `U`, along with the Index connecting `U` to `S`). - -The code in the `for` loop iterates over the diagonal elements of the `S` tensor (which are the singular values from the SVD), computes their squares to obtain the probabilities of observing the various states in the Schmidt basis (i.e. eigenvectors of the left-right bipartition reduced density matrices), and puts them into the von Neumann entanglement entropy formula ``S_\text{vN} = - \sum_{n} p_{n} \log{p_{n}}``. - -## Sampling from an MPS - -A matrix product state (MPS) can be viewed as defining a probability distribution -through the Born rule, as is the case when the MPS represents a quantum wavefunction. -To sample from the distribution defined by an MPS, you can use the function `sample` -provided in ITensor. For an MPS `psi` call to `sample(psi)` returns a random -sample from the distribution defined by `psi`. (Note that each sample is drawn anew -and not from a Markov chain seeded by a previous sample; this is possible because -the algorithm for sampling MPS is a `perfect' sampling algorithm with no autocorrelation.) - -In more detail, say we have a set of `N` site indices `s` and define a random MPS -with these sites: -```@example sample_mps; continued=true -using ITensors, ITensorMPS - -N = 10 # number of sites -d = 3 # dimension of each site -chi = 16 # bond dimension of the MPS -s = siteinds(d,N) -psi = random_mps(s;linkdims=chi) -``` - -We can now draw some samples from this MPS as - -```@example sample_mps -v1 = sample(psi) -v2 = sample(psi) -v3 = sample(psi) -println(v1) -println(v2) -println(v3) -``` - -The integers in each of the samples represent settings of each of the MPS indices -in the "computational basis". - -For reasons of efficiency, the `sample` function requires the MPS to be in orthogonal -form, orthogonalized to the first site. If it is not already in this form, it -can be brought into orthogonal form by calling `psi = orthogonalize(psi, 1)`. - - -## Write and Read an MPS or MPO to Disk with HDF5 - -!!! info - Make sure to install the HDF5 package to use this feature. (Run `julia> ] add HDF5` in the Julia REPL console.) - -**Writing an MPS to an HDF5 File** - -Let's say you have an MPS `psi` which you have made or obtained -from a calculation. To write it to an HDF5 file named "myfile.h5" -you can use the following pattern: - -```julia -using HDF5 -f = h5open("myfile.h5","w") -write(f,"psi",psi) -close(f) -``` - -Above, the string "psi" can actually be any string you want such as "MPS psi" -or "Result MPS" and doesn't have to have the same name as the reference `psi`. -Closing the file `f` is optional and you can also write other objects to the same -file before closing it. - -**Reading an MPS from an HDF5 File** - -Say you have an HDF5 file "myfile.h5" which contains an MPS stored as a dataset with the -name "psi". (Which would be the situation if you wrote it as in the example above.) -To read this ITensor back from the HDF5 file, use the following pattern: - -```julia -using HDF5 -f = h5open("myfile.h5","r") -psi = read(f,"psi",MPS) -close(f) -``` - -Many functions which involve MPS, such as the `dmrg` function or the `OpSum` system -require that you use an array of site indices which match the MPS. So when reading in -an MPS from disk, do not construct a new array of site indices. Instead, you can -obtain them like this: `sites = siteinds(psi)`. - -So for example, to create an MPO from an OpSum which has the same site indices -as your MPS `psi`, do the following: - -```julia -using ITensors, ITensorMPS - -os = OpSum() -# Then put operators into os... - -sites = siteinds(psi) # Get site indices from your MPS -H = MPO(os,sites) - -# Compute -energy_psi = inner(psi',H,psi) -``` - -Note the `MPS` argument to the read function, which tells Julia which read function -to call and how to interpret the data stored in the HDF5 dataset named "psi". In the -future we might lift the requirement of providing the type and have it be detected -automatically from the data stored in the file. - - -**Writing and Reading MPOs** - -To write or read MPOs to or from HDF5 files, just follow the examples above but use -the type `MPO` when reading an MPO from the file instead of the type `MPS`. - diff --git a/docs/src/examples/Physics.md b/docs/src/examples/Physics.md deleted file mode 100644 index 751e76e4f4..0000000000 --- a/docs/src/examples/Physics.md +++ /dev/null @@ -1,549 +0,0 @@ -# Physics (SiteType) System Examples - -## Obtaining a Predefined Operator - -Given an Index carrying a "physical" tag such as "Qubit", "S=1/2", "Boson", etc. -there are a set of pre-defined operators for each tag. The entire set of operators -can be found in the section [SiteTypes Included with ITensor](@ref). - -If you have an Index `s` carrying a "S=1/2" tag, for example, you can obtain the "Sz" -operator like this: -```julia -using ITensors, ITensorMPS - -op("Sz",s) -``` - -Usually indices with physical tags come from an array of indices returned from the `siteinds` function -```julia -using ITensors, ITensorMPS - -N = 10 -sites = siteinds("S=1/2",N) -``` -in which case one might want the "Sz" operator on site 4 -```julia -using ITensors, ITensorMPS -Sz4 = op("Sz",sites[4]) -``` - -## Make a Custom Operator from a Matrix - -The `op` function can be passed any matrix, as long as it has the correct dimensions, -and it will make this into an ITensor representing the operator with the corresponding -matrix elements. - -For example, if we have a two-dimensional Index `s` we could make the "Sz" operator ourselves from -the matrix -```julia -M = [1/2 0 ; 0 -1/2] -``` -by calling -```julia -using ITensors, ITensorMPS -Sz = op(M,s) -``` - - -## [Making a Custom op Definition](@id custom_op) - -The function `op` is used to obtain operators defined for a -given "site type". ITensor includes pre-defined site types such -as "S=1/2", "S=1", "Electron" and others. Or you can define your own site type -as discussed in detail in the code examples further below. - -**Extending op Function Definitions** - -Perhaps the most common part of the site type system one wishes to extend -are the various `op` or `op!` function overloads which allow code like - -```julia -using ITensors, ITensorMPS -s = siteind("S=1/2") -Sz = op("Sz",s) -``` - -to automatically create the ``S^z`` operator for an Index `s` based on the -`"S=1/2"` tag it carries. A major reason to define such `op` overloads -is to allow the OpSum system to recognize new operator names, as -discussed more below. - -Let's see how to introduce a new operator name into the ITensor site type -system for this existing site type of `"S=1/2"`. The operator we will -introduce is the projector onto the up spin state ``P_\uparrow`` which -we will denote with the string `"Pup"`. - -As a matrix acting on the space ``\{ |\!\uparrow\rangle, |\!\downarrow\rangle \}``, -the ``P_\uparrow`` operator is given by - -```math -\begin{aligned} - -P_\uparrow &= -\begin{bmatrix} - 1 & 0 \\ - 0 & 0 \\ -\end{bmatrix} - -\end{aligned} -``` - -To add this operator to the ITensor `op` system, we just need to introduce the following -code - -```julia -using ITensors, ITensorMPS -ITensors.op(::OpName"Pup",::SiteType"S=1/2") = - [1 0 - 0 0] -``` - -This code can be defined anywhere, such as in your own personal application code and does -not have to be put into the ITensor library source code. - -Note that we have to name the function `ITensors.op` and not just `op` so that it overloads -other functions of the name `op` inside the ITensors module. - -Having defined the above code, we can now do things like - -```julia -using ITensors, ITensorMPS -s = siteind("S=1/2") -Pup = op("Pup",s) -``` - -to obtain the `"Pup"` operator for our `"S=1/2"` Index `s`. Or we can do a similar -thing for an array of site indices: - -```julia -using ITensors, ITensorMPS -N = 40 -s = siteinds("S=1/2",N) -Pup1 = op("Pup",s[1]) -Pup3 = op("Pup",s[3]) -``` - -Note that for the `"Qudit"`/`"Boson"` site types, you have to define your overload -of `op` with the dimension of the local Hilbert space, for example: -```julia -using ITensors, ITensorMPS -function ITensors.op(::OpName"P1", ::SiteType"Boson", d::Int) - o = zeros(d, d) - o[1, 1] = 1 - return o -end -``` -Alternatively you could use Julia's [array comprehension](https://docs.julialang.org/en/v1/manual/arrays/#man-comprehensions) syntax: -```julia -using ITensors, ITensorMPS -ITensors.op(::OpName"P1", ::SiteType"Boson", d::Int) = - [(i == j == 1) ? 1.0 : 0.0 for i in 1:d, j in 1:d] -``` - -**Using Custom Operators in OpSum** - -A key use of these `op` system extensions is allowing additional operator names to -be recognized by the OpSum system for constructing matrix product operator (MPO) -tensor networks. With the code above defining the `"Pup"` operator, we are now -allowed to use this operator name in any OpSum code involving `"S=1/2"` site -indices. - -For example, we could now make an OpSum involving our custom operator such as: - -```julia -using ITensors, ITensorMPS -N = 100 -sites = siteinds("S=1/2",N) -os = OpSum() -for n=1:N - os += "Pup",n -end -P = MPO(os,sites) -``` - -This code makes an MPO `P` which is just the sum of a spin-up projection operator -acting on every site. - - -## Making a Custom state Definition - -The function `state` is used to define states (single-site wavefunctions) -that sites can be in. For example, the "Qubit" site type includes -definitions for the "0" and "1" states as well as the "+" (eigenstate of X operator) -state. The "S=1/2" site type includes definitions for the "Up" and "Dn" (down) states. - -Say we want to define a new state for the "Electron" site type called "+", which has -the meaning of one electron with its spin in the +X direction. First let's review -the existing state definitions: -```julia -using ITensors, ITensorMPS -ITensors.state(::StateName"Emp", ::SiteType"Electron") = [1.0, 0, 0, 0] -ITensors.state(::StateName"Up", ::SiteType"Electron") = [0.0, 1, 0, 0] -ITensors.state(::StateName"Dn", ::SiteType"Electron") = [0.0, 0, 1, 0] -ITensors.state(::StateName"UpDn", ::SiteType"Electron") = [0.0, 0, 0, 1] -``` -As we can see, the four settings of an "Electron" index correspond to the states -``|0\rangle, |\uparrow\rangle, |\downarrow\rangle, |\uparrow\downarrow\rangle``. - -So we can define our new state "+" as follows: -```julia -ITensors.state(::StateName"+", ::SiteType"Electron") = [0, 1/sqrt(2), 1/sqrt(2), 0] -``` -which makes the state -```math -|+\rangle = \frac{1}{\sqrt{2}} |\uparrow\rangle + \frac{1}{\sqrt{2}} |\downarrow\rangle -``` - -Having defined this overload of `state`, if we have an Index of type "Electron" -we can obtain our new state for it by doing -```julia -using ITensors, ITensorMPS -s = siteind("Electron") -plus = state("+",s) -``` -We can also use this new state definition in other ITensor features such as -the MPS constructor taking an array of state names. - - -## Make a Custom Local Hilbert Space / Physical Degree of Freedom - -ITensor provides support for a range of common local Hilbert space types, -or physical degrees of freedom, such as S=1/2 and S=1 spins; spinless and spinful -fermions; and more. - -However, there can be many cases where you need to make custom -degrees of freedom. You might be working with an -exotic system, such as ``Z_N`` parafermions for example, or need -to customize other defaults provided by ITensor. - -In ITensor, such a customization is done by overloading functions -on specially designated Index tags. -Below we give an brief introduction by example of how to make -such custom Index site types in ITensor. -Other code formulas following this one explain how to build on this -example to expand the capabilities of your custom site type such as -adding support for quantum number (QN) conservation and defining -custom mappings of strings to states. - -Throughout we will focus on the example of ``S=3/2`` spins. These -are spins taking the ``S^z`` values of ``+3/2,+1/2,-1/2,-3/2``. -So as tensor indices, they are indices of dimension 4. - -The key operators we will make for this example are ``S^z``, ``S^+``, -and ``S^-``, which are defined as: - -```math -\begin{aligned} -S^z &= -\begin{bmatrix} -3/2 & 0 & 0 & 0 \\ - 0 & 1/2 & 0 & 0 \\ - 0 & 0 &-1/2 & 0 \\ - 0 & 0 & 0 &-3/2\\ -\end{bmatrix} \\ - -S^+ & = -\begin{bmatrix} - 0 & \sqrt{3} & 0 & 0 \\ - 0 & 0 & 2 & 0 \\ - 0 & 0 & 0 & \sqrt{3} \\ - 0 & 0 & 0 & 0 \\ -\end{bmatrix} \\ - -S^- & = -\begin{bmatrix} - 0 & 0 & 0 & 0 \\ - \sqrt{3} & 0 & 0 & 0 \\ - 0 & 2 & 0 & 0 \\ - 0 & 0 & \sqrt{3} & 0 \\ -\end{bmatrix} \\ -\end{aligned} -``` - -**Code Preview** - -First let's see the minimal code needed to define and use this new -``S=3/2`` site type, then we will discuss what each part of -the code is doing. - -```julia -using ITensors, ITensorMPS - -ITensors.space(::SiteType"S=3/2") = 4 - -ITensors.op(::OpName"Sz",::SiteType"S=3/2") = - [+3/2 0 0 0 - 0 +1/2 0 0 - 0 0 -1/2 0 - 0 0 0 -3/2] - -ITensors.op(::OpName"S+",::SiteType"S=3/2") = - [0 √3 0 0 - 0 0 2 0 - 0 0 0 √3 - 0 0 0 0] - -ITensors.op(::OpName"S-",::SiteType"S=3/2") = - [0 0 0 0 - √3 0 0 0 - 0 2 0 0 - 0 0 √3 0] - -``` - -Now let's look at each part of the code above. - -**The SiteType** - -The most important aspect of this code is a special type, known as a `SiteType`, -which is a type made from a string. The string of interest here will be an Index -tag. In the code above, the `SiteType` we are using is - -```julia -SiteType"S=3/2" -``` - -What is the purpose of a `SiteType`? The answer is that we would like to be -able to select different functions to call on an ITensor Index based on what tags -it has, but that is not directly possible in Julia or indeed most languages. -However, if we can map a tag -to a type in the Julia type system, we can create function overloads for that type. -ITensor does this for certain functions for you, and we will discuss a few of these -functions below. So if the code encounters an Index such as `Index(4,"S=3/2")` it can -call these functions which are specialized for indices carrying the `"S=3/2"` tag. - -**The space Function** - -One of the overloadable `SiteType` functions is `space`, whose job is to -describe the vector space corresponding to that site type. For our -`SiteType"S=3/2"` overload of `space`, which gets called for any Index -carrying the `"S=3/2"` tag, the definition is - -```julia -using ITensors, ITensorMPS -ITensors.space(::SiteType"S=3/2") = 4 -``` - -Note that the function name is prepended with `ITensors.` before `space`. -This prefix makes sure the function is overloading other versions of the `space` -inside the `ITensors` module. - -The only information needed about the vector space of a `"S=3/2"` Index in -this example is that it is of dimension four. So the `space` function returns -the integer `4`. We will see in more advanced examples that the returned value -can instead be an array which specifies not only the dimension of a `"S=3/2"` -Index, but also additional subspace structure it has corresponding to quantum -numbers. - -After defining this `space` function, you can just write code like: - -```julia -using ITensors, ITensorMPS -s = siteind("S=3/2") -``` - -to obtain a single `"S=3/2"` Index, or write code like - -```julia -using ITensors, ITensorMPS -N = 100 -sites = siteinds("S=3/2",N) -``` - -to obtain an array of N `"S=3/2"` indices. The custom `space` function -will be used to determine the dimension of these indices, and the `siteind` -or `siteinds` functions provided by ITensor will help with extra things like -putting other Index tags that are conventional for site indices. - -**The op Function** - -The `op` function lets you define custom local operators associated -to the physical degrees of freedom of your `SiteType`. Then for example -you can use indices carrying your custom tag with OpSum and the -OpSum system will know how to automatically convert names of operators -such as `"Sz"` or `"S+"` into ITensors so that it can make an actual MPO. - -In our example above, we defined this function for the case of the `"Sz"` -operator as: - -```@example S32 -using ITensors, ITensorMPS -ITensors.op(::OpName"Sz",::SiteType"S=3/2") = - [+3/2 0 0 0 - 0 +1/2 0 0 - 0 0 -1/2 0 - 0 0 0 -3/2] -``` - -As you can see, the function is passed two objects: an `OpName` and a `SiteType`. -The strings `"Sz"` and `"S=3/2"` are also part of the type of these objects, and -have the meaning of which operator name we are defining and which site type these -operators are defined for. - -The body of this overload of `ITensors.op` constructs and returns a Julia matrix -which gives the matrix elements of the operator we are defining. - -Once this function is defined, and if you have an Index such as - -```@example S32; continued = true -s = Index(4,"S=3/2") -``` - -then, for example, you can get the `"Sz"` operator for this Index -and print it out by doing: - -```@example S32 -using ITensors, ITensorMPS -Sz = op("Sz",s) -println(Sz) -``` - -Again, through the magic of the `SiteType` -system, the ITensor library takes your Index, reads off its tags, -notices that one of them is `"S=3/2"`, and converts this into the type -`SiteType"S=3/2"` in order to call the specialized function `ITensors.op` defined above. - -You can use the `op` function yourself with a set of site indices created from -the `siteinds` function like this: - -```julia -using ITensors, ITensorMPS -N = 100 -sites = siteinds("S=3/2",N) -Sz1 = op("Sz",sites[1]) -Sp3 = op("S+",sites[3]) -``` - -Alternatively, you can write the lines of code above in the style -of `Sz1 = op("Sz",sites,1)`. - -This same `op` function is used inside of OpSum (formerly called AutoMPO) -when it converts its input into -an actual MPO. So by defining custom operator names you can pass any of these -operator names into OpSum and it will know how to use these operators. - -**Further Steps** - -See how the built-in site types are defined inside the ITensor library: -* [S=1/2 sites](https://github.com/ITensor/ITensors.jl/blob/main/src/lib/SiteTypes/src/sitetypes/spinhalf.jl) - Dimension 2 local Hilbert space. Similar to the `"Qubit"` site type, shares many of the same operator definitions. -* [Qubit sites](https://github.com/ITensor/ITensors.jl/blob/main/src/lib/SiteTypes/src/sitetypes/qubit.jl) - Dimension 2 local Hilbert space. Similar to the `"S=1/2"` site type, shares many of the same operator definitions. -* [S=1 sites](https://github.com/ITensor/ITensors.jl/blob/main/src/lib/SiteTypes/src/sitetypes/spinone.jl) - Dimension 3 local Hilbert space. -* [Fermion sites](https://github.com/ITensor/ITensors.jl/blob/main/src/lib/SiteTypes/src/sitetypes/fermion.jl) - Dimension 2 local Hilbert space. Spinless fermion site type. -* [Electron sites](https://github.com/ITensor/ITensors.jl/blob/main/src/lib/SiteTypes/src/sitetypes/electron.jl) - Dimension 4 local Hilbert space. Spinfull fermion site type. -* [tJ sites](https://github.com/ITensor/ITensors.jl/blob/main/src/lib/SiteTypes/src/sitetypes/tj.jl) - Dimension 3 local Hilbert space. Spinfull fermion site type but without a doubly occupied state in the Hilbert space. -* [Boson sites](https://github.com/ITensor/ITensors.jl/blob/main/src/lib/SiteTypes/src/sitetypes/boson.jl) - General d-dimensional local Hilbert space. Shares the same operator definitions as the `"Qudit"` site type. -* [Qudit sites](https://github.com/ITensor/ITensors.jl/blob/main/src/lib/SiteTypes/src/sitetypes/qudit.jl) - General d-dimensional local Hilbert space. Generalization of the `"Qubit"` site type, shares the same operator definitions as the ``Boson`` site type. - - -## Make a Custom Local Hilbert Space with QNs - -In the previous example above, we discussed the basic, -minimal code needed to define a custom local Hilbert space, using the example -of a ``S=3/2`` spin Hilbert space. In those examples, the `space` function -defining the vector space of a ``S=3/2`` spin only provides the dimension of -the space. But the Hilbert space of a ``S=3/2`` spin has additional structure, which -is that each of its four subspaces (each of dimension 1) can be labeled by -a different ``S^z`` quantum number. - -In this code formula we will include this extra quantum information in the -definition of the space of a ``S=3/2`` spin. - -**Code Preview** - -First let's see the minimal code needed to add the option for including -quantum numbers of our ``S=3/2`` site type, then we will discuss what each part of -the code is doing. - -```julia -using ITensors, ITensorMPS - -function ITensors.space(::SiteType"S=3/2"; - conserve_qns=false) - if conserve_qns - return [QN("Sz",3)=>1,QN("Sz",1)=>1, - QN("Sz",-1)=>1,QN("Sz",-3)=>1] - end - return 4 -end - -ITensors.op(::OpName"Sz",::SiteType"S=3/2") = - [+3/2 0 0 0 - 0 +1/2 0 0 - 0 0 -1/2 0 - 0 0 0 -3/2] - -ITensors.op(::OpName"S+",::SiteType"S=3/2") = - [0 √3 0 0 - 0 0 2 0 - 0 0 0 √3 - 0 0 0 0] - -ITensors.op(::OpName"S-",::SiteType"S=3/2") = - [0 0 0 0 - √3 0 0 0 - 0 2 0 0 - 0 0 √3 0] - - -``` - -Now let's look at each part of the code above. - -**The space function** - -In the previous code example above, we discussed -that the function `space` tells the ITensor library the basic information about how -to construct an Index associated with a special Index tag, in this case the tag `"S=3/2"`. -As in that code formula, if the user does not request that quantum numbers be included -(the case `conserve_qns=false`) then all that the `space` function returns is the number -4, indicating that a `"S=3/2"` Index should be of dimension 4. - -But if the `conserve_qns` keyword argument gets set to `true`, the `space` function we -defined above returns an array of `QN=>Int` pairs. (The notation `a=>b` in Julia constructs -a `Pair` object.) Each pair in the array denotes a subspace. -The `QN` part of each pair says what quantum number the subspace has, and the integer following -it indicates the dimension of the subspace. - -After defining the `space` function this way, you can write code like: - -```julia -using ITensors, ITensorMPS -s = siteind("S=3/2"; conserve_qns=true) -``` - -to obtain a single `"S=3/2"` Index which carries quantum number information. -The `siteind` function built into ITensor relies on your custom `space` function -to ask how to construct a `"S=3/2"` Index but also includes some other Index tags -which are conventional for all site indices. - -You can now also call code like: - -```julia -using ITensors, ITensorMPS -N = 100 -sites = siteinds("S=3/2",N; conserve_qns=true) -``` - -to obtain an array of N `"S=3/2"` indices which carry quantum numbers. - -**The op Function in the Quantum Number Case** - -Note that the `op` function overloads are exactly the same as for the -more basic case of defining an `"S=3/2"` Index type that does not carry -quantum numbers. There is no need to upgrade any of the `op` functions -for the QN-conserving case. -The reason is that all QN, block-sparse information -about an ITensor is deduced from the indices of the tensor, and setting elements -of such tensors does not require any other special code. - -However, only operators which have a well-defined QN flux---meaning they always -change the quantum number of a state they act on by a well-defined amount---can -be used in practice in the case of QN conservation. Attempting to build an operator, or any ITensor, -without a well-defined QN flux out of QN-conserving indices will result in a run time error. -An example of an operator that would lead to such an error would be the "Sx" spin operator -since it alternately increases ``S^z`` or decreases ``S^z`` depending on the state it acts -on, thus it does not have a well-defined QN flux. But it is perfectly fine to define an -`op` overload for the "Sx" operator and to make this operator when working with dense, -non-QN-conserving ITensors or when ``S^z`` is not conserved. - - diff --git a/docs/src/examples/mps_element.png b/docs/src/examples/mps_element.png deleted file mode 100644 index 648069ffd9..0000000000 Binary files a/docs/src/examples/mps_element.png and /dev/null differ diff --git a/docs/src/examples/mps_expect.png b/docs/src/examples/mps_expect.png deleted file mode 100644 index f36ecdb855..0000000000 Binary files a/docs/src/examples/mps_expect.png and /dev/null differ diff --git a/docs/src/examples/mps_from_tensor.png b/docs/src/examples/mps_from_tensor.png deleted file mode 100644 index 5298f8689b..0000000000 Binary files a/docs/src/examples/mps_from_tensor.png and /dev/null differ diff --git a/docs/src/examples/mps_onesite_figures/operator_app_mps.png b/docs/src/examples/mps_onesite_figures/operator_app_mps.png deleted file mode 100644 index bfd82589bb..0000000000 Binary files a/docs/src/examples/mps_onesite_figures/operator_app_mps.png and /dev/null differ diff --git a/docs/src/examples/mps_onesite_figures/operator_contract.png b/docs/src/examples/mps_onesite_figures/operator_contract.png deleted file mode 100644 index 45d222b5e5..0000000000 Binary files a/docs/src/examples/mps_onesite_figures/operator_contract.png and /dev/null differ diff --git a/docs/src/examples/mps_onesite_figures/updated_mps.png b/docs/src/examples/mps_onesite_figures/updated_mps.png deleted file mode 100644 index 35f33fe5ef..0000000000 Binary files a/docs/src/examples/mps_onesite_figures/updated_mps.png and /dev/null differ diff --git a/docs/src/examples/mps_zz_correlation.png b/docs/src/examples/mps_zz_correlation.png deleted file mode 100644 index ddc3ecd401..0000000000 Binary files a/docs/src/examples/mps_zz_correlation.png and /dev/null differ diff --git a/docs/src/examples/twosite_figures/gate_app_mps.png b/docs/src/examples/twosite_figures/gate_app_mps.png deleted file mode 100644 index 108b0987d8..0000000000 Binary files a/docs/src/examples/twosite_figures/gate_app_mps.png and /dev/null differ diff --git a/docs/src/examples/twosite_figures/gate_contract.png b/docs/src/examples/twosite_figures/gate_contract.png deleted file mode 100644 index 5242dc80c0..0000000000 Binary files a/docs/src/examples/twosite_figures/gate_contract.png and /dev/null differ diff --git a/docs/src/examples/twosite_figures/gate_gauge.png b/docs/src/examples/twosite_figures/gate_gauge.png deleted file mode 100644 index 6287c4ba95..0000000000 Binary files a/docs/src/examples/twosite_figures/gate_gauge.png and /dev/null differ diff --git a/docs/src/examples/twosite_figures/gate_svd.png b/docs/src/examples/twosite_figures/gate_svd.png deleted file mode 100644 index dc592a34ed..0000000000 Binary files a/docs/src/examples/twosite_figures/gate_svd.png and /dev/null differ diff --git a/docs/src/faq/DMRG.md b/docs/src/faq/DMRG.md deleted file mode 100644 index fa32ae4892..0000000000 --- a/docs/src/faq/DMRG.md +++ /dev/null @@ -1,242 +0,0 @@ -# Density Matrix Renormalization Group (DMRG) Frequently Asked Questions - -## Ensuring a DMRG calculation is converged - -While DMRG calculations can be extremely quick to converge in the best cases, -convergence can be slower for cases such as gapless systems or quasi-two-dimensional systems. -So it becomes important to know if a DMRG calculation is converged i.e. has been run -long enough with enough resources (large enough MPS bond dimension). - -Unfortunately **there is no automatic or bulletproof check for DMRG convergence**. -However, there are a number of reliable heuristics you can use to check convergence. -We list some of these with the most fundamental and important ones first: - -* Run your DMRG calculation on a **smaller system** and compare with another method, such - as an exact diagonalization. If the agreement is good, then gradually try larger - systems and see if the physical properties are roughly consistent and similar (i.e. - the density profile has similar features). - -* Make sure to check a **wide range of properties** - not just the energy. See if these - look plausible by plotting and visually inspecting them. For example: if your system has - left-right reflection symmetry, does the density or magnetization also have this symmetry? - If the ground state of your system is expected to have a total ``S^z`` of zero, does your - ground state have this property? - -* Make sure to run your DMRG calculation for **different numbers of sweeps** to see if - the results change. For example, if you run DMRG for 5 sweeps but are unsure of convergence, - try running it for 10 sweeps: is the energy the same or has it significantly decreased? - If 10 sweeps made a difference, try 20 sweeps. - -* Try setting the `eigsolve_krylovdim` keyword argument to a higher value (the default is 3). - This can be particularily helpful when the Hamiltonian is close to a critical point. - This may make slowly-converging calculations converge in fewer sweeps, but setting it - too high can make each sweep run slowly. - -* Inspect the the **DMRG output**. - The ITensor DMRG code reports the maximum bond or link dimension and maximum truncation error - after each sweep. (The maximums here mean over each DMRG substep making up one sweep.) - Is the maximum dimension or "maxlinkdim" reported by the DMRG output quickly reaching - and saturating the maxdim value you set for each sweep? Is the maximum truncation error - "maxerr" consistently reaching large values, larger than 1E-5? - Then it you may need to raise the maxdim parameter for your later sweeps, - so that DMRG is allowed to use a larger bond dimension and thus reach a better accuracy. - -* Compute the **energy variance** of an MPS to check whether it is an eigenstate. To do this - in ITensor, you can use the following code where `H` is your Hamiltonian MPO - and `psi` is the wavefunction you want to check: - - ```julia - H2 = inner(H,psi,H,psi) - E = inner(psi',H,psi) - var = H2-E^2 - @show var - ``` - - Here `var` is the quantity ``\langle H^2 \rangle - \langle H \rangle^2``. - The closer `var` is to zero, the more precisely `psi` is an eigenstate of `H`. Note - that this check does not ensure that `psi` is the ground state, but only one of the - eigenstates. - - -## Preventing DMRG from getting stuck in a local minimum - -While DMRG has very robust convergence properties when the initial MPS is close to the global -minimum, if it is far from the global minumum then there is _no guarantee_ that DMRG will -be able to find the true ground state. This problem is exacerbated for quantum number conserving -DMRG where the search space is more constrained. - -Thus it is very important to perform a number of checks to ensure that the result you -get from DMRG is actually converged. To learn about these checks, see the previous question. - -When DMRG is failing to converge, here are some of the steps you can take to improve things: - -* _The most important and useful technique_ is to turn on the **noise term** feature of DMRG. - To do this, just set the `noise` parameter of each sweep to a small, non-zero value, making - this value very small (1E-11, say) or zero by the last sweep. (Experiment with different - values on small systems to see which noise magnitudes help.) Here is an example of - defining DMRG accuracy or sweep parameters with a non-zero noise set for the first three sweeps: - - ```julia - nsweeps = 10 - maxdim = [100, 200, 400, 800, 1600] - cutoff = [1E-6] - noise = [1E-6, 1E-7, 1E-8, 0.0] - ... - energy, psi = dmrg(H,psi0; nsweeps, maxdim, cutoff, noise) - ``` - -* Try using a initial MPS with properties close to the ground state you are looking for. - For example, the ground state of a system of electrons typically has a density which is - spread out over the whole system. So if your initial state has all of the electrons bunched - up on the left-hand side only, it can take DMRG a very long time to converge. - -* Try using a random MPS with a modestly large bond dimension. ITensor offers a function - called [`random_mps`](@ref) which can be used to make random MPS in both the quantum number (QN) - conserving and non-QN conserving cases. Because random MPS have properties - which are "typical" of most ground states, they can be good initial states for DMRG. - -* Try DMRG on a closely related Hamiltonian for which convergence is easier to obtain - (be creative here: it could be your Hamiltonian with interactions turned off, or - with interactions only within, but not between, small local patches). Take the - output of this first calculation and use it as input for DMRG with the full Hamiltonian. - -* In stubborn cases, try other methods for finding the ground state which are slower, but - have a better chance of succeeding. A key example is imaginary time evolution, which - always reaches the ground state if (a) performed accurately on (b) a state which is - not orthogonal to the ground state. After doing some amount of imaginary time evolution, - use the resulting MPS as an initial state for DMRG obtain a higher-accuracy solution. - -## How to do periodic boundary condition DMRG - -The short answer to how to do fully periodic boundary condition DMRG in ITensor is that -you simply input a **periodic Hamiltonian** into our OpSum system and make the MPO -form of your Hamiltonian in the usual way. For example, for a chain of N sites with nearest-neighbor -interactions, you include a term that connects site 1 to site N. For a one-dimensional Ising model -chain Hamiltonian this would look like: - -``` -sites = siteinds("S=1/2",N) - -hterms = OpSum() -for j=1:(N-1) - hterms += "Sz",j,"Sz",j+1 -end -hterms += "Sz",1,"Sz",N # term 'wrapping' around the ring - -H = MPO(hterms,sites) -``` - -For two-dimensional DMRG calculations, where the most common approach is to use -periodic boundary conditions in the y-direction only, and not in the x-direction, -you do a similar step in making your OpSum input to ITensor DMRG: you include -terms wrapping around the periodic cylinder in the y direction but not in the x direction. - -However, fully periodic boundary conditions are only recommended for small systems -when absolutely needed, and in general are not recommended. For a longer discussion -of alternatives to using fully periodic boundaries, see the next section below. - -The reason fully periodic boundary conditions (periodic in x in 1D, and periodic in both x -and y in 2D) are not recommended in general is that the DMRG algorithm, as we are defining it -here, optimizes an **open-boundary MPS**. So if you input a periodic-boundary Hamiltonian, there -is a kind of "mismatch" that happens where you can still get the correct answer, but it -requires much more resources (a larger bond dimension and more sweeps) to get good -accuracy. There has been some research into "truly" periodic DMRG, [^Pippan] that is DMRG that -optimizes an MPS with a ring-like topology, but it is not widely used, is still an -open area of algorithm development, and is not currently available in ITensor. - - -## What boundary conditions should I choose: open, periodic, or infinite? - -One of the weaknesses of the density matrix renormalization group (DMRG), and its time-dependent or finite-temperature extensions, is that it works poorly with periodic boundary conditions. This stems from the fact that conventional DMRG optimizes over open-boundary matrix product state (MPS) wavefunctions whether or not the Hamiltonian includes periodic interactions. - -But this begs the question, when are periodic boundary conditions (PBC) really needed? DMRG offers -some compelling alternatives to PBC: - -* Use open boundary conditions (OBC). Though this introduces edge effects, the number of states needed - to reach a given accuracy is _significantly_ lower than with PBC (see next section below). - For gapped systems DMRG scales linearly with system size, meaning often one can study systems with many hundreds or even thousands of sites. Last but not least, open boundaries are often more natural. For studying systems which spontaneously break symmetry, adding "pinning" fields on the edge is often a very nice way to tip the balance toward a certain symmetry broken state while leaving the bulk unmodified. - -* Use smooth boundary conditions. The basic idea is to use OBC but - send the Hamiltonian parameters smoothly to zero at the boundary so that the system can not "feel" - the boundary. For certain systems this can significantly reduce edge effects.[^Smooth1][^Smooth2][^Smooth3] - -[^Smooth1]: [Smooth boundary conditions for quantum lattice systems](http://dx.doi.org/10.1103/PhysRevLett.71.4283), M. Vekic and Steven R. White, _Phys. Rev. Lett._ **71**, [4283](http://dx.doi.org/10.1103/PhysRevLett.71.4283) (1993) cond-mat/[9310053](http://arxiv.org/abs/cond-mat/9310053) - -[^Smooth2]: [Hubbard model with smooth boundary conditions](http://dx.doi.org/10.1103/PhysRevB.53.14552), M. Vekic and Steven R. White, _Phys. Rev. B_ **53**, [14552](http://dx.doi.org/10.1103/PhysRevB.53.14552) (1996) cond-mat/[9601009](http://arxiv.org/abs/cond-mat/9601009) - -[^Smooth3]: [Grand canonical finite-size numerical approaches: A route to measuring bulk properties in an applied field](http://link.aps.org/doi/10.1103/PhysRevB.86.041108), Chisa Hotta and Naokazu Shibata, _Phys. Rev. B_ **86**, [041108](http://link.aps.org/doi/10.1103/PhysRevB.86.041108) (2012) - - - -* Use "infinite boundary conditions", that is, use infinite DMRG in the form of an algorithm like iDMRG or VUMPS. This has a cost that can be even less than with OBC yet is completely free of finite-size effects. - -However, there are a handful of cases where PBC remains preferable despite the extra overhead. A few such cases are: - -* Benchmarking DMRG against another code that uses PBC, such as a Monte Carlo or exact diagonalization code. - -* Extracting the central charge of a critical one-dimensional system described by a CFT. In practice, using PBC can give an accurate central charge even for quite small systems by fitting the subsystem entanglement entropy to the CFT scaling form. - -* Checking for the presence or absence of topological effects. These could be edge effects (the Haldane - phase has a four-fold ground state degeneracy with OBC, but not with PBC), or could be related to some global topological sector that is ill-defined with PBC (e.g. periodic vs. antiperiodic boundary conditions for the transverse field Ising model). - -(Note that in the remaining discussion, by PBC I mean *fully periodic* boundary conditions in all directions. -For the case of DMRG applied to quasi-two-dimensional systems, it remains a good practice to use -periodic boundaries in the shorter direction, while still using open (or infinite) boundaries -in the longer direction along the DMRG/MPS path.) - -Below I discuss more about the problems with using PBC, as well as some misconceptions about when PBC seems necessary even though there are better alternatives. - -#### Drawbacks of Periodic Boundary Conditions - -Periodic boundary conditions are straightforward to implement in conventional DMRG. The simplest approach is to include a "long bond" directly connecting site 1 to site N in the Hamiltonian. However this -naive approach has a major drawback: if open-boundary DMRG achieves a given accuracy when keeping ``m`` states (bond dimension of size ``m``), then to reach the same accuracy with PBC one must keep closer to ``m^2`` states! The reason is that now every bond of the MPS not only carries local entanglement as with OBC, but also the entanglement between the first and last sites. (There is an alternative DMRG algorithm[^Pippan] for periodic systems which may have better scaling than the above approach but has not been widely applied and tested, as far as I am aware, especially for - 2D or critical systems .) - -[^Pippan]: [Efficient matrix-product state method for periodic boundary conditions](http://link.aps.org/doi/10.1103/PhysRevB.81.081103), P. Pippan, Steven R. White, and H.G. Evertz, _Phys. Rev. B_ **81**, [081103](http://link.aps.org/doi/10.1103/PhysRevB.81.081103) - -The change in scaling from ``m`` to ``m^2`` is a severe problem. -For example, many gapped one-dimensional systems only require about ``m=100`` to reach good accuracy -(truncation errors of less than 1E-9 or so). To reach the same accuracy with naive PBC would then -require using 10,000 states, which can easily fill the RAM of a typical desktop computer for a large enough system, not to mention the extra time needed to work with larger matrices. - -But poor scaling is not the only drawback of PBC. Systems that exhibit spontaneous symmetry breaking -are simple to work with under OBC, where one has the additional freedom of applying edge pinning terms -to drive the bulk into a specific symmetry sector. Using edge pinning reduces the bulk entanglement and makes measuring order parameters straightforward. Similarly one can use infinite DMRG to directly observe symmetry breaking effects. - -But under PBC, order parameters remain equal to zero and can only be accessed through correlation functions. Though using correlation functions is often presented as the "standard" or "correct" approach, such reasoning pre-supposes that PBC is the best choice. Recent work in the quantum Monte Carlo community demonstrates that open boundaries with pinning fields can actually be a superior approach.[^Assaad] - -[^Assaad]: [Pinning the Order: The Nature of Quantum Criticality in the Hubbard Model on Honeycomb Lattice](http://dx.doi.org/10.1103/PhysRevX.3.031010), Fakher F. Assaad and Igor F. Herbut, _Phys. Rev. X_ **3**, [031010](http://dx.doi.org/10.1103/PhysRevX.3.031010) - - -#### Cases Where Periodic BC Seems Necessary, But Open/Infinite BC Can be Better - -Below are some cases where periodic boundary conditions seem to be necessary at a first glance. -But in many of these cases, not only can open or infinite boundaries be just as successful, they -can even be the better choice. - -* _Measuring asymptotic properties of correlation functions_: much of our understanding of gapless one-dimensional systems comes from field-theoretic approaches which make specific predictions about asymptotic decays of various correlators. To test these predictions numerically, one must work with large, translationally invariant systems with minimal edge effects. Using fully periodic boundary conditions satisfies these criteria. However, a superior choice is to use infinite DMRG, which combines the much better scaling of open-boundary DMRG with the ability to measure correlators at _arbitrarily long_ distances by repeating the unit cell of the MPS wavefunction. Although truncating to a finite number of states imposes an effective correlation length on the system, this correlation length can reach many thousands of sites for quite moderate MPS bond dimensions. Karrasch and Moore took advantage of this fact to convincingly check the predictions of Luttinger liquid theory for one-dimensional systems of gapless fermions.[^Karrasch] - -[^Karrasch]: [Luttinger liquid physics from the infinite-system density matrix renormalization group](http://dx.doi.org/10.1103/PhysRevB.86.155156), C. Karrasch and J.E. Moore, _Phys. Rev. B_ **86**, [155156](http://dx.doi.org/10.1103/PhysRevB.86.155156) - -* _Studying two-dimensional topological order_: a hallmark of intrinsic topological order is the presence of a robust ground state degeneracy when the system is put on a torus. Also many topological phases have gapless edge states which can cause problems for numerical calculations. Thus one might think that fully periodic BC are the best choice for studying topological phases. However, topological phases have the same ground-state degeneracy on an infinite cylinder as they do on a torus.[^Zhang]. Cincio and Vidal exploited this fact to use infinite DMRG to study a variety of topological phases [^Cincio]. One part of their calculation did actually require obtaining ground states on a torus, but they accomplished this by taking a finite segment of an infinite MPS and connecting its ends. This approach does not give the true ground state of the torus but was sufficient for their calculation and was arguably closer to the true two-dimensional physics. - -[^Zhang]: [Quasiparticle statistics and braiding from ground-state entanglement](http://dx.doi.org/10.1103/PhysRevB.85.235151), Yi Zhang, Tarun Grover, Ari Turner, Masaki Oshkawa, and Ashvin Vishwanath, _Phys. Rev. B_ **85**, [235151](http://dx.doi.org/10.1103/PhysRevB.85.235151) - -[^Cincio]: [Characterizing Topological Order by Studying the Ground States on an Infinite Cylinder](http://link.aps.org/doi/10.1103/PhysRevLett.110.067208), L. Cincio and G. Vidal, _Phys. Rev. Lett._ **110**, [067208](http://link.aps.org/doi/10.1103/PhysRevLett.110.067208) - -* _Obtaining bulk gaps_: - DMRG has the ability to "target" low-lying excited states or to obtain such - states by constraining them to be orthogonal to the ground state. However, with OBC, - localized excitations can get stuck to the edges and not reveal the true bulk gap behavior. - Thus one may conclude that PBC is necessary. But using open or infinite boundaries remains - the better choice because they allow much higher accuracy. - - To deal with the presence of edges in OBC, one can use "restricted sweeping". Here one sweeps across the full system to obtain the ground state. Then, to obtain the first excited state one only sweeps through the full system to obtain the ground state. Then, to obtain the first excited state one only sweeps through the near the edges. This traps the particle in a "soft box" which still lets its wavefunction mix with the basis that describes the ground state outside the restricted sweeping region. - - Within infinite DMRG, boundary effects are rigorously absent if the calculation has converged. To compute bulk gaps one again uses a type of restricted sweeping known in the literature as "infinite boundary conditions". For more see the work by Phien, Vidal, and McCulloch.[^Phien] - -[^Phien]: [Infinite boundary conditions for matrix product state calculations](http://link.aps.org/doi/10.1103/PhysRevB.86.245107), Ho N. Phien, G. Vidal, and Ian P. McCulloch _Phys. Rev. B_ **86**, [245107](http://link.aps.org/doi/10.1103/PhysRevB.86.245107) - - -In conclusion, consider carefully whether you really need to use periodic boundary conditions, as they impose a steep computational cost within DMRG. Periodic BC can actually be worse for the very types of measurements where they are often presented as the best or "standard" choice. Many of the issues periodic boundaries circumvent can be avoided more elegantly by using infinite DMRG, or when that is not applicable, by using open boundary conditions with sufficient care. - diff --git a/docs/src/faq/QN.md b/docs/src/faq/QN.md deleted file mode 100644 index 7d66f22261..0000000000 --- a/docs/src/faq/QN.md +++ /dev/null @@ -1,24 +0,0 @@ -# Quantum Number Frequently Asked Questions - -## Can I mix different types of quantum numbers within the same system? - -Yes, you can freely mix quantum numbers (QNs) of different types. For example, -you can make the sites of your systems alternate between sites carrying -spin "Sz" QNs and fermion sites carrying particle number "Nf" QNs. The QNs will -not mix with each other and will separately be conserved to the original value -you set for your initial wavefunction. - -## How can I separately conserve QNs which have the same name? - -If you have two physically distinct types of sites, such as "Qudit" sites, but -which carry identically named QNs called "Number", and you want the qudit number -to be separately conserved within each type of site, -you must make the QN names different for the two types of sites. - -For example, the following line of code will make an array of site indices with the qudit number QN having the name "Number\_odd" on odd sites and "Number\_even" on even sites: -``` -sites = [isodd(n) ? siteind("Qudit", n; dim=10, conserve_qns=true, qnname_number="Number_odd") - : siteind("Qudit", n; dim=2, conserve_qns=true, qnname_number="Number_even") - for n in 1:2*L] -``` -(You may have to collapse the above code into a single line for it to run properly.) diff --git a/docs/src/faq/RelationshipToOtherLibraries.md b/docs/src/faq/RelationshipToOtherLibraries.md deleted file mode 100644 index 620c3a647e..0000000000 --- a/docs/src/faq/RelationshipToOtherLibraries.md +++ /dev/null @@ -1,13 +0,0 @@ -# Relationship of ITensor to other tensor libraries - -Here we will describe the relationship of ITensor to more traditional Julia Arrays or deep learning libraries like TensorFlow and PyTorch. There are a few things that distinguish ITensor from those approaches: - -1. ITensors have dimensions with labels that get passed around, which makes it simple to perform certain operations like contraction, addition, and tensor decompositions with a high level interface, independent of memory layout. This is along the same lines as Julia packages like [NamedDims.jl](https://github.com/invenia/NamedDims.jl) and [AxisArrays.jl](https://github.com/JuliaArrays/AxisArrays.jl) and libraries in Python like [xarray](https://xarray.pydata.org/en/stable/index.html), however I would argue that the ITensor approach is a little more sophisticated (the dimensions have more metadata which makes them easier to manipulate for different situations, random ids to help avoid name clashes, etc.). This design was inspired by the needs of tensor network algorithms, where there are many tensor dimensions in the computation (of which many of them are dynamically created during the calculation), but would be helpful for writing other algorithms too. - -2. The ITensor type has a dynamic high level interface, where the type itself is mutable and the data can be swapped out. This allows for conveniently allocating the data of an ITensor on the fly "as needed", which makes for a nicer, more flexible interface (like initializing an empty ITensor before a loop, and filling it with the correct data type when the first value is set), at the expense of a small overhead for accessing data in the ITensor. We have found this tradeoff is worth it, since we expect ITensors to be used for medium to large scale calculations where operations on the tensors like contraction, addition, and tensor decomposition dominate the cost of the calculation, and code can be designed with function barriers to speed up operations when data is being accessed repeatedly. - -3. Another feature that ITensor has that goes beyond what is available in standard Julia, TensorFlow, and PyTorch is tensors which are [symmetric under a group action](https://arxiv.org/pdf/1008.4774.pdf). The physical interpretation of these tensors are ones that have a conserved quantity (like a quantum state with a conserved number of particles), so that feature is more physics-oriented, but could have applications in other areas like machine learning as well. In practice, these tensors are block sparse, and have extra metadata on the dimensions labeling representations of the group. - -4. Based on the features above, the ITensor library provides high level implementations of tensor network algorithms (algebraic operations of very high dimensional tensors, such as addition, multiplication, and finding dominant eigenvectors). In general these algorithms can (and have been) written on top of other libraries like standard Julia Arrays/AD, PyTorch, or TensorFlow, but they might have various downsides (a less convenient interface for dealing with tensor operations, no support for the types of symmetric tensors we often need, limited support for tensors with complex numbers in the case of libraries like PyTorch, though perhaps that has improved since I last checked, etc.). - -Although ITensor has primarily focused on quantum physics and quantum computing applications, there is work using ITensor for machine learning applications (so far focused on applications of tensor networks to machine learning, so no neural network calculations yet as far as I know). In general, these different libraries (ITensor, Flux, PyTorch, TensorFlow) are biased towards their specific methods and application areas that they are used for the most: ITensor is more biased towards tensor network calculations and quantum physics/quantum computing applications, based on the available features and interface, while PyTorch and TensorFlow are more biased towards neural network calculations. However, our goal would be to provide more features to ITensor that would make it useful for neural network applications as well, such as better support for slicing operations. diff --git a/docs/src/getting_started/NextSteps.md b/docs/src/getting_started/NextSteps.md index a2eeaffb22..1e93fa07e1 100644 --- a/docs/src/getting_started/NextSteps.md +++ b/docs/src/getting_started/NextSteps.md @@ -1,12 +1,10 @@ # Next Steps -* Try one of the [Tutorials](@ref dmrg_tutorial) in the next section of the ITensor documentation - * Browse the [Code Examples](@ref itensor_examples). * Read the [ITensor Paper](https://www.scipost.org/SciPostPhysCodeb.4) for a long-form introduction to the design and main features of the ITensor library -* Read the [Advanced ITensor Usage Guide](@ref advanced_usage_guide) +* Refer to the [ITensorMPS.jl documentation](https://itensor.github.io/ITensorMPS.jl) for information on running finite MPS/MPO calculations such as DMRG. * More Julia language tutorials and resources - [From zero to Julia!](https://techytok.com/from-zero-to-julia/) diff --git a/docs/src/index.md b/docs/src/index.md index 8ca81e0b9f..0009cdd4dc 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -5,30 +5,33 @@ tensor network algorithms. | **Documentation**|**Citation**| |:----------------:|:----------:| -|[![docs](https://img.shields.io/badge/docs-latest-blue.svg)](https://itensor.github.io/ITensors.jl/dev/)|[![SciPost](https://img.shields.io/badge/SciPost-10.21468-blue.svg)](https://scipost.org/SciPostPhysCodeb.4) [![arXiv](https://img.shields.io/badge/arXiv-2007.14822-b31b1b.svg)](https://arxiv.org/abs/2007.14822)| +|[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://itensor.github.io/ITensors.jl/stable/) [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://itensor.github.io/ITensors.jl/dev/)|[![SciPost](https://img.shields.io/badge/SciPost-10.21468-blue.svg)](https://scipost.org/SciPostPhysCodeb.4) [![arXiv](https://img.shields.io/badge/arXiv-2007.14822-b31b1b.svg)](https://arxiv.org/abs/2007.14822)| |**Version**|**Download Statistics**|**Style Guide**|**License**| |:---------:|:---------------------:|:-------------:|:---------:| |[![version](https://juliahub.com/docs/ITensors/version.svg)](https://juliahub.com/ui/Packages/ITensors/P3pqL)|[![ITensor Downloads](https://img.shields.io/badge/dynamic/json?url=http%3A%2F%2Fjuliapkgstats.com%2Fapi%2Fv1%2Fmonthly_downloads%2FITensors&query=total_requests&suffix=%2Fmonth&label=Downloads)](http://juliapkgstats.com/pkg/ITensors)|[![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle)|[![license](https://img.shields.io/badge/License-Apache_2.0-blue.svg)](https://github.com/ITensor/ITensors.jl/blob/main/LICENSE)| -The source code for ITensor can be found [on Github](https://github.com/ITensor/ITensors.jl). +The source code for ITensors.jl can be found [on Github](https://github.com/ITensor/ITensors.jl). -Additional documentation can be found on the ITensor website [itensor.org](https://itensor.org/). +Additional documentation can be found on the ITensor website [itensor.org](https://itensor.org). An ITensor is a tensor whose interface is independent of its memory layout. ITensor indices are objects which carry extra information and which 'recognize' each other (compare equal to each other). -The ITensor library also includes composable and extensible -algorithms for optimizing and transforming tensor networks, such as -matrix product state and matrix product operators, such as -the DMRG algorithm. +The [ITensorMPS.jl library](https://github.com/ITensor/ITensorMPS.jl) +includes composable and extensible algorithms for optimizing and transforming +tensor networks, such as matrix product state and matrix product operators, such as +the DMRG algorithm. If you are looking for information on running finite MPS/MPO +calculations such as DMRG, take a look at the [ITensorMPS.jl documentation](https://itensor.github.io/ITensorMPS.jl). Development of ITensor is supported by the Flatiron Institute, a division of the Simons Foundation. ## News +- March 22, 2025: As part of the latest release of ITensors.jl (v0.8.3), all documentation related to MPS/MPO functionality has been moved to the [ITensorMPS.jl documentation](https://itensor.github.io/ITensorMPS.jl). + - February 22, 2025: Please note that there were issues installing the latest version of ITensors.jl (ITensors.jl v0.8) in older versions of Julia v1.10 and v1.11 ([https://github.com/ITensor/ITensors.jl/issues/1618](https://github.com/ITensor/ITensors.jl/issues/1618), [https://itensor.discourse.group/t/typeparameteraccessors-not-found-error-on-julia-v-1-10-0/2260](https://itensor.discourse.group/t/typeparameteraccessors-not-found-error-on-julia-v-1-10-0/2260)). This issue has been fixed in [NDTensors.jl v0.4.4](https://github.com/ITensor/ITensors.jl/pull/1623), so please try updating your packages if you are using older versions of Julia v1.10 or v1.11 and running into issues installing ITensors.jl. - February 3, 2025: ITensors.jl v0.8 has been released. This release should not be breaking to the average user using documented features of the library. This removes internal submodules that held experimental code for rewriting the internals of NDTensors.jl/ITensors.jl, which have now been turned into separate packages for future development. It is marked as breaking since ITensorMPS.jl was making use of some of that experimental code, and will be updated accordingly. Also note that it fixes an issue that existed in some more recent versions of NDTensors.jl v0.3/ITensors.jl v0.7 where loading ITensors.jl in combination with some packages like LinearMaps.jl caused very long load/compile times ([https://itensor.discourse.group/t/linearmaps-and-itensors-incompatibility/2216](https://itensor.discourse.group/t/linearmaps-and-itensors-incompatibility/2216)), so if you are seeing that issue when using ITensors.jl v0.7 you should upgrade to this version. @@ -271,61 +274,3 @@ hastags(s, "Site") = true hastags(s, "n=1") = true i1 == i = false ``` - -### DMRG Calculation - -DMRG is an iterative algorithm for finding the dominant -eigenvector of an exponentially large, Hermitian matrix. -It originates in physics with the purpose of finding -eigenvectors of Hamiltonian (energy) matrices which model -the behavior of quantum systems. - -```julia -using ITensors, ITensorMPS -let - # Create 100 spin-one indices - N = 100 - sites = siteinds("S=1",N) - - # Input operator terms which define - # a Hamiltonian matrix, and convert - # these terms to an MPO tensor network - # (here we make the 1D Heisenberg model) - os = OpSum() - for j=1:N-1 - os += "Sz",j,"Sz",j+1 - os += 0.5,"S+",j,"S-",j+1 - os += 0.5,"S-",j,"S+",j+1 - end - H = MPO(os,sites) - - # Create an initial random matrix product state - psi0 = random_mps(sites) - - # Plan to do 5 passes or 'sweeps' of DMRG, - # setting maximum MPS internal dimensions - # for each sweep and maximum truncation cutoff - # used when adapting internal dimensions: - nsweeps = 5 - maxdim = [10,20,100,100,200] - cutoff = 1E-10 - - # Run the DMRG algorithm, returning energy - # (dominant eigenvalue) and optimized MPS - energy, psi = dmrg(H,psi0; nsweeps, maxdim, cutoff) - println("Final energy = $energy") - - nothing -end - -# output - -After sweep 1 energy=-137.954199761732 maxlinkdim=9 maxerr=2.43E-16 time=9.356 -After sweep 2 energy=-138.935058943878 maxlinkdim=20 maxerr=4.97E-06 time=0.671 -After sweep 3 energy=-138.940080155429 maxlinkdim=92 maxerr=1.00E-10 time=4.522 -After sweep 4 energy=-138.940086009318 maxlinkdim=100 maxerr=1.05E-10 time=11.644 -After sweep 5 energy=-138.940086058840 maxlinkdim=96 maxerr=1.00E-10 time=12.771 -Final energy = -138.94008605883985 -``` -You can find more examples of running `dmrg` and related algorithms [here](https://github.com/ITensor/ITensorMPS.jl/tree/main/examples). - diff --git a/docs/src/tutorials/DMRG.md b/docs/src/tutorials/DMRG.md deleted file mode 100644 index 9ff95b823d..0000000000 --- a/docs/src/tutorials/DMRG.md +++ /dev/null @@ -1,127 +0,0 @@ -# [DMRG Tutorial](@id dmrg_tutorial) - -The [density matrix renormalization group (DMRG)](https://tensornetwork.org/mps/algorithms/dmrg/) -is an algorithm for computing eigenstates -of Hamiltonians (or extremal eigenvectors of large, Hermitian matrices). -It computes these eigenstates in the -[matrix product state (MPS)](https://tensornetwork.org/mps/) format. - -Let's see how to set up and run a DMRG calculation using the ITensor library. -We will be interested in finding the ground state of the quantum Hamiltonian -``H`` given by: - -```math -H = \sum_{j=1}^{N-1} \mathbf{S}_{j} \cdot \mathbf{S}_{j+1} = \sum_{j=1}^{N-1} S^z_{j} S^z_{j+1} + \frac{1}{2} S^+_{j} S^-_{j+1} + \frac{1}{2} S^-_{j} S^+_{j+1} -``` - -This Hamiltonian is known as the one-dimensional Heisenberg model and we will -take the spins to be ``S=1`` spins (spin-one spins). We will consider -the case of ``N=100`` and plan to do five sweeps of DMRG (five passes over the system). - -**ITensor DMRG Code** - -Let's look at an entire, working ITensor code that will do this calculation then -discuss the main steps. If you need help running the code below, see the getting -started page on [Running ITensor and Julia Codes](@ref). - -```julia -using ITensors, ITensorMPS -let - N = 100 - sites = siteinds("S=1",N) - - os = OpSum() - for j=1:N-1 - os += "Sz",j,"Sz",j+1 - os += 1/2,"S+",j,"S-",j+1 - os += 1/2,"S-",j,"S+",j+1 - end - H = MPO(os,sites) - - psi0 = random_mps(sites;linkdims=10) - - nsweeps = 5 - maxdim = [10,20,100,100,200] - cutoff = [1E-10] - - energy,psi = dmrg(H,psi0;nsweeps,maxdim,cutoff) - - return -end -``` - - -**Steps of The Code** - -The first two lines - -```@example siteinds; continued=true -using ITensors, ITensorMPS -N = 100 -sites = siteinds("S=1",N) -``` - -tells the function `siteinds` to make an array of ITensor [Index](https://itensor.github.io/ITensors.jl/stable/IndexType.html) objects which -have the properties of ``S=1`` spins. This means their dimension will be 3 and -they will carry the `"S=1"` tag, which will enable the next part of the code to know -how to make appropriate operators for them. - -Try printing out some of these indices to verify their properties: - -```@example siteinds -@show sites[1] -``` - -The next part of the code builds the Hamiltonian: - -```julia -os = OpSum() -for j=1:N-1 - os += "Sz",j,"Sz",j+1 - os += 1/2,"S+",j,"S-",j+1 - os += 1/2,"S-",j,"S+",j+1 -end -H = MPO(os,sites) -``` - -An `OpSum` is an object which accumulates Hamiltonian terms such as `"Sz",1,"Sz",2` -so that they can be summed afterward into a matrix product operator (MPO) tensor network. -The line of code `H = MPO(os,sites)` constructs the Hamiltonian in the MPO format, with -physical indices given by the array `sites`. - -The line - -```julia -psi0 = random_mps(sites;linkdims=10) -``` - -constructs an MPS `psi0` which has the physical indices `sites` and a bond dimension of 10. -It is made by a random quantum circuit that is reshaped into an MPS, so that it will have as generic and unbiased properties as an MPS of that size can have. -This choice can help prevent the DMRG calculation from getting stuck in a local minimum. - -The lines - -```julia -nsweeps = 5 -maxdim = [10,20,100,100,200] -cutoff = [1E-10] -``` - -define the number of DMRG sweeps (five) we will instruct the code to do, as well as the -parameters that will control the speed and accuracy of the DMRG algorithm within -each sweep. The array `maxdim` limits the maximum MPS bond dimension allowed during -each sweep and `cutoff` defines the truncation error goal of each sweep (if fewer values are -specified than sweeps, the last value is used for all remaining sweeps). - -Finally the call - -```julia -energy,psi = dmrg(H,psi0;nsweeps,maxdim,cutoff) -``` - -runs the DMRG algorithm included in ITensor, using `psi0` as an -initial guess for the ground state wavefunction. The optimized MPS `psi` and -its eigenvalue `energy` are returned. - -After the `dmrg` function returns, you can take the returned MPS `psi` and do further calculations with it, such as measuring local operators or computing entanglement entropy. - diff --git a/docs/src/tutorials/MPSTimeEvolution.md b/docs/src/tutorials/MPSTimeEvolution.md deleted file mode 100644 index 58886a6c6b..0000000000 --- a/docs/src/tutorials/MPSTimeEvolution.md +++ /dev/null @@ -1,189 +0,0 @@ -# MPS Time Evolution - -An important application of [matrix product state (MPS)](https://tensornetwork.org/mps/) -tensor networks in physics is computing the time evolution of a quantum state under the dynamics -of a Hamiltonian ``H``. An accurate, efficient, and simple way to time evolve a matrix product state (MPS) is by using a Trotter decomposition of the time evolution operator ``U(t) = e^{-i H t}``. - -The technique we will use is "time evolving block decimation" (TEBD). -More simply it is just the idea of decomposing the time-evolution operator into a circuit of -quantum 'gates' (two-site unitaries) using the Trotter-Suzuki approximation and applying these gates in -a controlled way to an MPS. - -Let's see how to set up and run a TEBD calculation using ITensor. - -The Hamiltonian ``H`` we will use is the one-dimensional Heisenberg model -which is given by: - -```math -\begin{aligned} -H & = \sum_{j=1}^{N-1} \mathbf{S}_{j} \cdot \mathbf{S}_{j+1} \\ -& = \sum_{j=1}^{N-1} S^z_{j} S^z_{j+1} + \frac{1}{2} S^+_{j} S^-_{j+1} + \frac{1}{2} S^-_{j} S^+_{j+1} -\end{aligned} -``` - -**The TEBD Method** - -When the Hamiltonian, like the one above, is a sum of local terms, - -```math -H = \sum_j h_{j,j+1} -``` - -where ``h_{j,j+1}`` acts on sites j and j+1, -then a Trotter decomposition that is particularly well suited for use -with MPS techniques is - -```math -e^{-i \tau H} \approx e^{-i h_{1,2} \tau/2} e^{-i h_{2,3} \tau/2} \cdots e^{-i h_{N-1,N} \tau/2} -e^{-i h_{N-1,N} \tau/2} e^{-i h_{N-2,N-1} \tau/2} \cdots e^{-i h_{1,2} \tau/2} + O(\tau^3) -``` - -Note the factors of two in each exponential. Each factored exponential is known as a -Trotter "gate". - -We can visualize the resulting circuit that will be applied to the MPS as follows: - -![](trotter_tevol.png) - -The error in the above decomposition is of order ``\tau^3``, so that will be the error -accumulated _per time step_. Because of the time-step error, one takes ``\tau`` to be -small and then applies the above set of operators to an MPS as a single sweep, then -does a number ``(t/\tau)`` of sweeps to evolve for a total time ``t``. The total error -will therefore scale as ``\tau^2`` with this scheme, though other sources of error may -dominate for long times, or very small ``\tau``, such as truncation errors. - -Let's take a look at the code to apply these Trotter gates to an MPS to -time evolve it. Then we will break down the steps of the code in more detail. - - -**ITensor TEBD Time Evolution Code** - -Let's look at an entire, working ITensor code that will do this calculation then -discuss the main steps. (If you need help running the code below, see the getting -started page on running ITensor codes.) - -```julia -using ITensors, ITensorMPS - -let - N = 100 - cutoff = 1E-8 - tau = 0.1 - ttotal = 5.0 - - # Make an array of 'site' indices - s = siteinds("S=1/2", N; conserve_qns=true) - - # Make gates (1,2),(2,3),(3,4),... - 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(-im * tau / 2 * hj) - push!(gates, Gj) - end - # Include gates in reverse order too - # (N,N-1),(N-1,N-2),... - append!(gates, reverse(gates)) - - # Initialize psi to be a product state (alternating up and down) - psi = MPS(s, n -> isodd(n) ? "Up" : "Dn") - - c = div(N, 2) # center site - - # Compute and print at each time step - # then apply the gates to go to the next time - for t in 0.0:tau:ttotal - Sz = expect(psi, "Sz"; sites=c) - println("$t $Sz") - - t≈ttotal && break - - psi = apply(gates, psi; cutoff) - normalize!(psi) - end - - return -end -``` - -**Steps of The Code** - -First we setsome parameters, like the system size N and time step ``\tau`` to use. - -The line `s = siteinds("S=1/2",N;conserve_qns=true)` defines an array of -spin 1/2 tensor indices (Index objects) which will be the site or physical -indices of the MPS. - -Next we make an empty array `gates = ITensor[]` that will hold ITensors -that will be our Trotter gates. Inside the `for n=1:N-1` loop that follows -the lines - -```julia -hj = op("Sz",s1) * op("Sz",s2) + - 1/2 * op("S+",s1) * op("S-",s2) + - 1/2 * op("S-",s1) * op("S+",s2) -``` - -call the `op` function which reads the "S=1/2" tag on our site indices -(sites j and j+1) and which then knows that we want the spin 1/ -2 version of the "Sz", "S+", and "S-" operators. -The `op` function returns these operators as ITensors and we -tensor product and add them together to compute the operator ``h_{j,j+1}`` -defined as - -```math -h_{j,j+1} = S^z_j S^z_{j+1} + \frac{1}{2} S^+_j S^-_{j+1} + \frac{1}{2} S^-_j S^+_{j+1} -``` - -which we call `hj` in the code. - -To make the corresponding Trotter gate `Gj` we exponentiate `hj` times -a factor ``-i \tau/2`` and then append or push this onto the end of the -gate array `gates`. - -```julia -Gj = exp(-im * tau/2 * hj) -push!(gates,Gj) -``` - -Having made the gates for bonds (1,2),(2,3),(3,4), etc. we still need -to append the gates in reverse order to complete the correct Trotter -formula. Here we can conveniently do that by just calling the Julia -`append!` function and supply a reversed version of the array of -gates we have made so far. This can -be done in a single line of code `append!(gates,reverse(gates))`. - -The line of code `psi = MPS(s, n -> isodd(n) ? "Up" : "Dn")` -initializes our MPS `psi` as a product state of alternating -up and down spins. - -To carry out the time evolution we loop over -the range of times from 0.0 to `ttotal` in steps of `tau`, -using the Julia range notation `0.0:tau:ttotal` to easily -set up this loop as `for t in 0.0:tau:ttotal`. - -Inside the loop, we use the `expect` function to measure -the expected value of the `"Sz"` operator on the center -site. - -To evolve the MPS to the next time, we call the function - -```julia -psi = apply(gates, psi; cutoff) -``` - -which applies the array of ITensors called `gates` to our current -MPS `psi`, truncating the MPS at each step using the truncation -error threshold supplied as the variable `cutoff`. - -The `apply` function is smart enough to determine which site indices -each gate has, and then figure out where to apply it to our -MPS. It automatically handles truncating the MPS and can -even handle non-nearest-neighbor gates, though that -feature is not used in this example. - diff --git a/docs/src/tutorials/QN_DMRG.md b/docs/src/tutorials/QN_DMRG.md deleted file mode 100644 index 3d3f593372..0000000000 --- a/docs/src/tutorials/QN_DMRG.md +++ /dev/null @@ -1,182 +0,0 @@ -# Quantum Number Conserving DMRG - -An important technique in DMRG calculations of quantum Hamiltonians -is the conservation of _quantum numbers_. Examples of these are the -total number of particles of a model of fermions, or the total of all -``S^z`` components of a system of spins. Not only can conserving quantum -numbers make DMRG calculations run more quickly and use less memory, but -it can be important for simulating physical systems with conservation -laws and for obtaining ground states in different symmetry sectors. -Note that ITensor currently only supports Abelian quantum numbers. - -#### Necessary Changes - -Setting up a quantum-number conserving DMRG calculation in ITensor requires -only very small changes to a DMRG code. The main changes are: - -1. using tensor indices (`Index` objects) which carry quantum number (QN) information to build your Hamiltonian and initial state -2. initializing your MPS to have well-defined total quantum numbers - -Importantly, _the total QN of your state throughout the calculation will -remain the same as the initial state passed to DMRG_. -The total QN of your state is not set separately, but determined -implicitly from the initial QN of the state when it is first constructed. - -Of course, your Hamiltonian should conserve all of the QN's that you would -like to use. If it doesn't, you will get an error when you try to construct -it out of the QN-enabled tensor indices. - -#### Making the Changes - -Let's see how to make these two changes to the -[DMRG Tutorial](@ref dmrg_tutorial) code from the previous section. -At the end, we will put together these changes for a complete, working code. - -**Change 1: QN Site Indices** - -To make change (1), we will change the line - -```julia -sites = siteinds("S=1",N) -``` - -by setting the `conserve_qns` keyword argument to `true`: - -```julia -sites = siteinds("S=1",N; conserve_qns=true) -``` - -Setting `conserve_qns=true` tells the `siteinds` function to conserve -every possible quantum number associated to the site -type (which is `"S=1"` in this example). For ``S=1`` spins, this will turn on -total-``S^z`` conservation. -(For other site types that conserve multiple QNs, there are specific keyword -arguments available to track just a subset of conservable QNs.) -We can check this by printing out some of the site indices, and seeing that the -subspaces of each `Index` are labeled by QN values: - -```julia -@show sites[1] -@show sites[2] -``` - -Sample output: - -``` - sites[1] = (dim=3|id=794|"S=1,Site,n=1") - 1: QN("Sz",2) => 1 - 2: QN("Sz",0) => 1 - 3: QN("Sz",-2) => 1 - sites[2] = (dim=3|id=806|"S=1,Site,n=2") - 1: QN("Sz",2) => 1 - 2: QN("Sz",0) => 1 - 3: QN("Sz",-2) => 1 -``` - -In the sample output above, note that in addition to the dimension of these indices being 3, each of the three settings of the Index have a unique QN associated to them. The number after the QN on each line is the dimension of that subspace, which is 1 for each subspace of the Index objects above. Note also that `"Sz"` quantum numbers in ITensor are measured in units of ``1/2``, so `QN("Sz",2)` corresponds to ``S^z=1`` in conventional physics units. - -**Change 2: Initial State** - -To make change (2), instead of constructing the initial MPS `psi0` to be an arbitrary, random MPS, we will make it a specific state with a well-defined total ``S^z``. -So we will replace the line - -```julia -psi0 = random_mps(sites;linkdims=10) -``` - -by the lines - -```julia -state = [isodd(n) ? "Up" : "Dn" for n=1:N] -psi0 = MPS(sites,state) -``` - -The first line of the new code above makes an array of strings which -alternate between `"Up"` and `"Dn"` on odd and even numbered sites. -These names `"Up"` and `"Dn"` are special values associated to the `"S=1"` -site type which indicate up and down spin values. The second line takes -the array of site Index objects `sites` and the array of strings `state` -and returns an MPS which is a product state (classical, unentangled state) -with each site's state given by the strings in the `state` array. -In this example, `psi0` will be a Neel state with alternating up and down -spins, so it will have a total ``S^z`` of zero. We could check this by -computing the quantum-number flux of `psi0` - -```julia -@show flux(psi0) -# Output: flux(psi0) = QN("Sz",0) -``` - -!!! info "Setting Other Total QN Values" - - The above example shows the case of setting a total "Sz" quantum - number of zero, since the initial state alternates between "Up" - and "Dn" on every site with an even number of sites. - - To obtain other total QN values, just set the initial state to - be one which has the total QN you want. To be concrete - let's take the example of a system with `N=10` sites of - ``S=1`` spins. - - For example if you want a total "Sz" of +20 (= `QN("Sz",20)`) in ITensor units, - or ``S^z=10`` in physical units, for a system with 10 sites, - use the initial state: - ```julia - state = ["Up" for n=1:N] - psi0 = MPS(sites,state) - ``` - Or to initialize this 10-site system to have a total "Sz" of +16 - in ITensor units (``S^z=8`` in physical units): - ```julia - state = ["Dn","Up","Up","Up","Up","Up","Up","Up","Up","Up"] - psi0 = MPS(sites,state) - ``` - would work (as would any `state` with one "Dn" and nine "Up"'s - in any order). - Or you could initialize to a total "Sz" of +18 - in ITensor units (``S^z=9`` in physical units) as - ```julia - state = ["Z0","Up","Up","Up","Up","Up","Up","Up","Up","Up"] - psi0 = MPS(sites,state) - ``` - where "Z0" refers to the ``S^z=0`` state of a spin-one spin. - - Finally, the same kind of logic as above applies to other - physical site types, whether "S=1/2", "Electron", - etc. - -#### Putting it All Together - -Let's take the [DMRG Tutorial](@ref dmrg_tutorial) code -from the previous section and make the changes discussed above, -to turn it into a code which conserves the total ``S^z`` quantum -number throughout the DMRG calculation. The resulting code is: - -```julia -using ITensors, ITensorMPS -let - N = 100 - sites = siteinds("S=1",N;conserve_qns=true) - - os = OpSum() - for j=1:N-1 - os += "Sz",j,"Sz",j+1 - os += 1/2,"S+",j,"S-",j+1 - os += 1/2,"S-",j,"S+",j+1 - end - H = MPO(os,sites) - - state = [isodd(n) ? "Up" : "Dn" for n=1:N] - psi0 = MPS(sites,state) - @show flux(psi0) - - nsweeps = 5 - maxdim = [10,20,100,100,200] - cutoff = [1E-10] - - energy, psi = dmrg(H,psi0; nsweeps, maxdim, cutoff) - - return -end -``` - diff --git a/docs/src/tutorials/tebd.jl b/docs/src/tutorials/tebd.jl deleted file mode 100644 index 7dfb452b7d..0000000000 --- a/docs/src/tutorials/tebd.jl +++ /dev/null @@ -1,46 +0,0 @@ -using ITensors, ITensorMPS - -let - N = 100 - cutoff = 1E-8 - tau = 0.1 - ttotal = 5.0 - - # Make an array of 'site' indices - s = siteinds("S=1/2", N; conserve_qns=true) - - # Make gates (1,2),(2,3),(3,4),... - 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 - # Include gates in reverse order too - # (N,N-1),(N-1,N-2),... - append!(gates, reverse(gates)) - - # Initialize psi to be a product state (alternating up and down) - psi = MPS(s, n -> isodd(n) ? "Up" : "Dn") - - c = div(N, 2) # center site - - # Compute and print at each time step - # then apply the gates to go to the next time - for t in 0.0:tau:ttotal - Sz = expect(psi, "Sz"; sites=c) - println("$t $Sz") - - t ≈ ttotal && break - - psi = apply(gates, psi; cutoff) - normalize!(psi) - end - - return nothing -end diff --git a/docs/src/tutorials/trotter_tevol.png b/docs/src/tutorials/trotter_tevol.png deleted file mode 100644 index bcf10ce93a..0000000000 Binary files a/docs/src/tutorials/trotter_tevol.png and /dev/null differ