diff --git a/Project.toml b/Project.toml index e064de7..8139685 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensorMPS" uuid = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" authors = ["Matthew Fishman ", "Miles Stoudenmire "] -version = "0.3.9" +version = "0.3.10" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/README.md b/README.md index a7f4b37..888fefa 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,13 @@ # ITensorMPS.jl -[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://ITensor.github.io/ITensorMPS.jl/stable/) -[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://ITensor.github.io/ITensorMPS.jl/dev/) +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://itensor.github.io/ITensorMPS.jl/stable/) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://itensor.github.io/ITensorMPS.jl/dev/) [![Build Status](https://github.com/ITensor/ITensorMPS.jl/actions/workflows/Tests.yml/badge.svg?branch=main)](https://github.com/ITensor/ITensorMPS.jl/actions/workflows/Tests.yml?query=branch%3Amain) [![Coverage](https://codecov.io/gh/ITensor/ITensorMPS.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/ITensor/ITensorMPS.jl) [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) [![Aqua](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) -Finite MPS and MPO methods based on the Julia version of [ITensor](https://www.itensor.org) ([ITensors.jl](https://github.com/ITensor/ITensors.jl)). See the [ITensors.jl documentation](https://itensor.github.io/ITensors.jl/dev/) for more details. +Finite MPS and MPO methods based on the Julia version of [ITensor](https://www.itensor.org) ([ITensors.jl](https://github.com/ITensor/ITensors.jl)). See the [ITensorMPS.jl documentation](https://itensor.github.io/ITensorMPS.jl) for more details. ## News @@ -31,6 +31,64 @@ This release introduces a new (experimental) function `expand` for performing gl ITensorMPS.jl v0.2 has been released, which is a breaking release. It updates to using ITensorTDVP.jl v0.4, which has a number of breaking changes to the `tdvp`, `linsolve`, and `dmrg_x` functions. See the [ITensorTDVP.jl v0.4 release notes](https://github.com/ITensor/ITensorTDVP.jl/blob/main/README.md#itensortdvpjl-v04-release-notes) for details. +## Example 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 in 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). + --- *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* diff --git a/docs/Project.toml b/docs/Project.toml index cc02678..b50e658 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,11 @@ [deps] -ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" +ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" + +[compat] +Documenter = "1.9.0" +ITensorMPS = "0.3.9" +ITensors = "0.8.1" +Literate = "2.20.1" diff --git a/docs/make.jl b/docs/make.jl index 50b0182..cc0e4ee 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,18 +1,51 @@ -using ITensorMPS: ITensorMPS +using ITensorMPS +using ITensors using Documenter: Documenter, DocMeta, deploydocs, makedocs DocMeta.setdocmeta!(ITensorMPS, :DocTestSetup, :(using ITensorMPS); recursive=true) +DocMeta.setdocmeta!(ITensors, :DocTestSetup, :(using ITensors); recursive=true) include("make_index.jl") makedocs(; - modules=[ITensorMPS], + # Allows using ITensors.jl docstrings in ITensorMPS.jl documentation: + # https://github.com/JuliaDocs/Documenter.jl/issues/1734 + modules=[ITensorMPS, ITensors], authors="ITensor developers and contributors", sitename="ITensorMPS.jl", format=Documenter.HTML(; canonical="https://ITensor.github.io/ITensorMPS.jl", edit_link="main", assets=String[] ), - pages=["Home" => "index.md"], + pages=[ + "Home" => "index.md", + "Tutorials" => [ + "DMRG" => "tutorials/DMRG.md", + "Quantum Number Conserving DMRG" => "tutorials/QN_DMRG.md", + "MPS Time Evolution" => "tutorials/MPSTimeEvolution.md", + ], + "Code Examples" => [ + "MPS and MPO Examples" => "examples/MPSandMPO.md", + "DMRG Examples" => "examples/DMRG.md", + "Physics (SiteType) System Examples" => "examples/Physics.md", + ], + "Documentation" => [ + "MPS and MPO" => "MPSandMPO.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", + ], + "Frequently Asked Questions" => + ["DMRG FAQs" => "faq/DMRG.md", "Quantum Number (QN) FAQs" => "faq/QN.md"], + "HDF5 File Formats" => "HDF5FileFormats.md", + ], warnonly=true, ) diff --git a/docs/src/DMRG.md b/docs/src/DMRG.md new file mode 100644 index 0000000..cc15808 --- /dev/null +++ b/docs/src/DMRG.md @@ -0,0 +1,5 @@ +# DMRG + +```@docs +dmrg +``` diff --git a/docs/src/DMRGObserver.md b/docs/src/DMRGObserver.md new file mode 100644 index 0000000..8a56cfc --- /dev/null +++ b/docs/src/DMRGObserver.md @@ -0,0 +1,43 @@ +# 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/HDF5FileFormats.md b/docs/src/HDF5FileFormats.md new file mode 100644 index 0000000..c9c329f --- /dev/null +++ b/docs/src/HDF5FileFormats.md @@ -0,0 +1,30 @@ +## [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 new file mode 100644 index 0000000..f8396c8 --- /dev/null +++ b/docs/src/IncludedSiteTypes.md @@ -0,0 +1,388 @@ +# 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/MPSandMPO.md b/docs/src/MPSandMPO.md new file mode 100644 index 0000000..5474bbc --- /dev/null +++ b/docs/src/MPSandMPO.md @@ -0,0 +1,162 @@ +# 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/Observer.md b/docs/src/Observer.md new file mode 100644 index 0000000..66affcf --- /dev/null +++ b/docs/src/Observer.md @@ -0,0 +1,200 @@ +# [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 new file mode 100644 index 0000000..89234e2 --- /dev/null +++ b/docs/src/OpSum.md @@ -0,0 +1,14 @@ +# OpSum + +## Description + +```@docs +OpSum +``` + +## Methods + +```@docs +add! +MPO(::OpSum,::Vector{<:Index}) +``` diff --git a/docs/src/ProjMPO.md b/docs/src/ProjMPO.md new file mode 100644 index 0000000..8010775 --- /dev/null +++ b/docs/src/ProjMPO.md @@ -0,0 +1,23 @@ +# 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 new file mode 100644 index 0000000..8a713b2 --- /dev/null +++ b/docs/src/ProjMPOSum.md @@ -0,0 +1,22 @@ +# 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/SiteType.md b/docs/src/SiteType.md new file mode 100644 index 0000000..42f4d11 --- /dev/null +++ b/docs/src/SiteType.md @@ -0,0 +1,17 @@ +# 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 new file mode 100644 index 0000000..889caa7 --- /dev/null +++ b/docs/src/Sweeps.md @@ -0,0 +1,26 @@ +# 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 new file mode 100644 index 0000000..48d51fe --- /dev/null +++ b/docs/src/examples/DMRG.md @@ -0,0 +1,552 @@ +# 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 new file mode 100644 index 0000000..c8ee559 --- /dev/null +++ b/docs/src/examples/MPSandMPO.md @@ -0,0 +1,454 @@ +# 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 new file mode 100644 index 0000000..751e76e --- /dev/null +++ b/docs/src/examples/Physics.md @@ -0,0 +1,549 @@ +# 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 new file mode 100644 index 0000000..648069f Binary files /dev/null and b/docs/src/examples/mps_element.png differ diff --git a/docs/src/examples/mps_expect.png b/docs/src/examples/mps_expect.png new file mode 100644 index 0000000..f36ecdb Binary files /dev/null and b/docs/src/examples/mps_expect.png differ diff --git a/docs/src/examples/mps_from_tensor.png b/docs/src/examples/mps_from_tensor.png new file mode 100644 index 0000000..5298f86 Binary files /dev/null and b/docs/src/examples/mps_from_tensor.png 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 new file mode 100644 index 0000000..bfd8258 Binary files /dev/null and b/docs/src/examples/mps_onesite_figures/operator_app_mps.png differ diff --git a/docs/src/examples/mps_onesite_figures/operator_contract.png b/docs/src/examples/mps_onesite_figures/operator_contract.png new file mode 100644 index 0000000..45d222b Binary files /dev/null and b/docs/src/examples/mps_onesite_figures/operator_contract.png differ diff --git a/docs/src/examples/mps_onesite_figures/updated_mps.png b/docs/src/examples/mps_onesite_figures/updated_mps.png new file mode 100644 index 0000000..35f33fe Binary files /dev/null and b/docs/src/examples/mps_onesite_figures/updated_mps.png differ diff --git a/docs/src/examples/mps_zz_correlation.png b/docs/src/examples/mps_zz_correlation.png new file mode 100644 index 0000000..ddc3ecd Binary files /dev/null and b/docs/src/examples/mps_zz_correlation.png differ diff --git a/docs/src/examples/twosite_figures/gate_app_mps.png b/docs/src/examples/twosite_figures/gate_app_mps.png new file mode 100644 index 0000000..108b098 Binary files /dev/null and b/docs/src/examples/twosite_figures/gate_app_mps.png differ diff --git a/docs/src/examples/twosite_figures/gate_contract.png b/docs/src/examples/twosite_figures/gate_contract.png new file mode 100644 index 0000000..5242dc8 Binary files /dev/null and b/docs/src/examples/twosite_figures/gate_contract.png differ diff --git a/docs/src/examples/twosite_figures/gate_gauge.png b/docs/src/examples/twosite_figures/gate_gauge.png new file mode 100644 index 0000000..6287c4b Binary files /dev/null and b/docs/src/examples/twosite_figures/gate_gauge.png differ diff --git a/docs/src/examples/twosite_figures/gate_svd.png b/docs/src/examples/twosite_figures/gate_svd.png new file mode 100644 index 0000000..dc592a3 Binary files /dev/null and b/docs/src/examples/twosite_figures/gate_svd.png differ diff --git a/docs/src/faq/DMRG.md b/docs/src/faq/DMRG.md new file mode 100644 index 0000000..fa32ae4 --- /dev/null +++ b/docs/src/faq/DMRG.md @@ -0,0 +1,242 @@ +# 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 new file mode 100644 index 0000000..7d66f22 --- /dev/null +++ b/docs/src/faq/QN.md @@ -0,0 +1,24 @@ +# 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/index.md b/docs/src/index.md index 61cd7ea..57224ab 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -4,14 +4,14 @@ EditURL = "../../examples/README.jl" # ITensorMPS.jl -[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://ITensor.github.io/ITensorMPS.jl/stable/) -[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://ITensor.github.io/ITensorMPS.jl/dev/) +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://itensor.github.io/ITensorMPS.jl/stable/) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://itensor.github.io/ITensorMPS.jl/dev/) [![Build Status](https://github.com/ITensor/ITensorMPS.jl/actions/workflows/Tests.yml/badge.svg?branch=main)](https://github.com/ITensor/ITensorMPS.jl/actions/workflows/Tests.yml?query=branch%3Amain) [![Coverage](https://codecov.io/gh/ITensor/ITensorMPS.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/ITensor/ITensorMPS.jl) [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) [![Aqua](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) -Finite MPS and MPO methods based on the Julia version of [ITensor](https://www.itensor.org) ([ITensors.jl](https://github.com/ITensor/ITensors.jl)). See the [ITensors.jl documentation](https://itensor.github.io/ITensors.jl/dev/) for more details. +Finite MPS and MPO methods based on the Julia version of [ITensor](https://www.itensor.org) ([ITensors.jl](https://github.com/ITensor/ITensors.jl)). See the [ITensorMPS.jl documentation](https://itensor.github.io/ITensorMPS.jl) for more details. ## News @@ -35,6 +35,64 @@ This release introduces a new (experimental) function `expand` for performing gl ITensorMPS.jl v0.2 has been released, which is a breaking release. It updates to using ITensorTDVP.jl v0.4, which has a number of breaking changes to the `tdvp`, `linsolve`, and `dmrg_x` functions. See the [ITensorTDVP.jl v0.4 release notes](https://github.com/ITensor/ITensorTDVP.jl/blob/main/README.md#itensortdvpjl-v04-release-notes) for details. +## Example 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. + +````@example index +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 in 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). + --- *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* diff --git a/docs/src/tutorials/DMRG.md b/docs/src/tutorials/DMRG.md new file mode 100644 index 0000000..9ff95b8 --- /dev/null +++ b/docs/src/tutorials/DMRG.md @@ -0,0 +1,127 @@ +# [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 new file mode 100644 index 0000000..58886a6 --- /dev/null +++ b/docs/src/tutorials/MPSTimeEvolution.md @@ -0,0 +1,189 @@ +# 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 new file mode 100644 index 0000000..3d3f593 --- /dev/null +++ b/docs/src/tutorials/QN_DMRG.md @@ -0,0 +1,182 @@ +# 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 new file mode 100644 index 0000000..7dfb452 --- /dev/null +++ b/docs/src/tutorials/tebd.jl @@ -0,0 +1,46 @@ +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 new file mode 100644 index 0000000..bcf10ce Binary files /dev/null and b/docs/src/tutorials/trotter_tevol.png differ diff --git a/examples/Project.toml b/examples/Project.toml index e513705..8e962c6 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -1,2 +1,5 @@ [deps] ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" + +[compat] +ITensorMPS = "0.3" diff --git a/examples/README.jl b/examples/README.jl index ca0a530..d439fa0 100644 --- a/examples/README.jl +++ b/examples/README.jl @@ -1,13 +1,13 @@ # # ITensorMPS.jl # -# [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://ITensor.github.io/ITensorMPS.jl/stable/) -# [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://ITensor.github.io/ITensorMPS.jl/dev/) +# [![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://itensor.github.io/ITensorMPS.jl/stable/) +# [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://itensor.github.io/ITensorMPS.jl/dev/) # [![Build Status](https://github.com/ITensor/ITensorMPS.jl/actions/workflows/Tests.yml/badge.svg?branch=main)](https://github.com/ITensor/ITensorMPS.jl/actions/workflows/Tests.yml?query=branch%3Amain) # [![Coverage](https://codecov.io/gh/ITensor/ITensorMPS.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/ITensor/ITensorMPS.jl) # [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) # [![Aqua](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) -# Finite MPS and MPO methods based on the Julia version of [ITensor](https://www.itensor.org) ([ITensors.jl](https://github.com/ITensor/ITensors.jl)). See the [ITensors.jl documentation](https://itensor.github.io/ITensors.jl/dev/) for more details. +# Finite MPS and MPO methods based on the Julia version of [ITensor](https://www.itensor.org) ([ITensors.jl](https://github.com/ITensor/ITensors.jl)). See the [ITensorMPS.jl documentation](https://itensor.github.io/ITensorMPS.jl) for more details. # # ## News # @@ -30,3 +30,59 @@ # #### Breaking changes # # ITensorMPS.jl v0.2 has been released, which is a breaking release. It updates to using ITensorTDVP.jl v0.4, which has a number of breaking changes to the `tdvp`, `linsolve`, and `dmrg_x` functions. See the [ITensorTDVP.jl v0.4 release notes](https://github.com/ITensor/ITensorTDVP.jl/blob/main/README.md#itensortdvpjl-v04-release-notes) for details. +# +# ## Example 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. + +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 in 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).