From 637dd2409f27ede41aa916822ea8acb4cd557a9e Mon Sep 17 00:00:00 2001 From: Ben Corbett <32752943+corbett5@users.noreply.github.com> Date: Mon, 26 Feb 2024 21:00:36 -0700 Subject: [PATCH] Corbett/types and docs (#9) * Specify type of MPO and pass OpCacheVec as strings * Docs and examples --- README.md | 135 ++++++++++++--------------- examples/electronic-structure.jl | 152 ++++++++++++++++++++++++------- examples/fermi-hubbard.jl | 132 ++++++++++++++++++++++----- examples/haldane-shastry.jl | 37 ++++++++ examples/intro.jl | 36 ++++++++ src/MPOConstruction.jl | 50 +++++++--- src/OpIDSum.jl | 31 +++++-- 7 files changed, 423 insertions(+), 150 deletions(-) create mode 100644 examples/haldane-shastry.jl create mode 100644 examples/intro.jl diff --git a/README.md b/README.md index 863d399..66681de 100644 --- a/README.md +++ b/README.md @@ -6,7 +6,7 @@ [![Coverage](https://codecov.io/gh/ITensor/ITensorMPOConstruction.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/ITensor/ITensorMPOConstruction.jl) [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) -A fast algorithm for constructing a Matrix Product Operator (MPO) from a sum of local operators. This is a replacement for `MPO(os::OpSum, sites::Vector{<:Index})`. In all cases examined so far this algorithm constructs an MPO with a smaller (or equal) bond dimension faster than the competition. +A fast algorithm for constructing a Matrix Product Operator (MPO) from a sum of local operators. This is a replacement for `MPO(os::OpSum, sites::Vector{<:Index})`. In all cases examined so far this algorithm constructs an MPO with a smaller (or equal) bond dimension faster than the competition. All runtimes below are taken from a single sample on a 2021 MacBook Pro with the M1 Max CPU and 32GB of memory. ## Installation @@ -17,15 +17,15 @@ julia> using Pkg; Pkg.add(url="https://github.com/ITensor/ITensorMPOConstruction ## Constraints -This algorithm shares same constraints as ITensor's default algorithm. +This algorithm shares the same constraints as ITensors' default algorithm. -* The operator must be expressed as a sum of products of single site operators. For example a CNOT could not appear in the sum since it is a two site operator. -* When dealing with Fermionic systems the parity of each term in the sum must be even. That is the combined number of creation and annihilation operators must be even. +1. The operator must be expressed as a sum of products of single site operators. For example a CNOT could not appear in the sum since it is a two site operator. +2. When dealing with Fermionic systems the parity of each term in the sum must be even. That is the combined number of creation and annihilation operators in each term in the sum must be even. There are also two additional constraints: -* Each term in the sum of products representation can only have a single operator acting on a site. For example a term such as $\mathbf{X}^{(1)} \mathbf{X}^{(1)}$ is not allowed. However, there is a pre-processing utility that can automatically replace $\mathbf{X}^{(1)} \mathbf{X}^{(1)}$ with $\mathbf{I}^{(1)}$. This is not a hard requirement for the algorithm, just a simplification to improve performance. Many operators of interest can be easily expressed in a form where only a single operator acts on each site in a term. -* When constructing a quantum number (QN) conserving operator the total flux of the operator must be zero. It would be trivial to remove this constraint. +3. Each term in the sum of products representation can only have a single operator acting on a site. For example a term such as $\mathbf{X}^{(1)} \mathbf{X}^{(1)}$ is not allowed. However, there is a pre-processing utility that can automatically replace $\mathbf{X}^{(1)} \mathbf{X}^{(1)}$ with $\mathbf{I}^{(1)}$. This is not a hard requirement for the algorithm but a simplification to improve performance. +4. When constructing a quantum number conserving operator the total flux of the operator must be zero. It would be trivial to remove this constraint. ## `MPO_new` @@ -36,8 +36,22 @@ function MPO_new(os::OpSum, sites::Vector{<:Index}; kwargs...)::MPO ``` The optional keyword arguments are -* `tol::Real`: The tolerance used in the sparse QR decomposition. The default value is calculated separately for each QR decomposition, it is almost always a good value. -* `basisOpCacheVec::OpCacheVec`: A list of operators to use as a basis for each site. The operators on at each site are expressed as one of these basis operators. +* `basisOpCacheVec`: A list of operators to use as a basis for each site. The operators on each site are expressed as one of these basis operators. +* `tol::Real`: The tolerance used in the sparse QR decomposition (which is done by SPQR). The default value is the SPQR default which is calculated separately for each QR decomposition. If you want a MPO that is accurate up to floating point errors the default tolerance should work well. If instead you want to compress the MPO the value `tol` will differ from the `cutoff` passed to `ITensor.MPO` since the truncation method is completely different. If you want to replicate the same truncation behavior first construct the MPO with a suitably small (or default) `tol` and then use `ITensors.truncate!`. + +## Examples: Fermi-Hubbard Hamiltonian in Real Space + +The one dimensional Fermi-Hubbard Hamiltonian with periodic boundary conditions on $N$ sites can be expressed in real space as + +$$ +\mathcal{H} = -t \sum_{i = 1}^N \sum_{\sigma \in (\uparrow, \downarrow)} \left( c^\dagger_{i, \sigma} c_{i + 1, \sigma} + c^\dagger_{i, \sigma} c_{i - 1, \sigma} \right) + U \sum_{i = 1}^N n_{i, \uparrow} n_{i, \downarrow} \ , +$$ + +where the periodic boundary conditions enforce that $c_k = c_{k + N}$. For this Hamiltonian all that needs to be done to switch over to using ITensorMPOConstruction is switch `MPO(os, sites)` to `MPO_New(os, sites)`. + +https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/fermi-hubbard.jl#L4-L24 + +For $N = 1000$ both ITensors and ITensorMPOConstruction can construct an MPO of bond dimension 10 in under two seconds. For a compelling reason to use ITensorMPOConstruction we need to look at a more complicated Hamiltonian. ## Examples: Fermi-Hubbard Hamiltonian in Momentum Space @@ -47,86 +61,57 @@ $$ \mathcal{H} = \sum_{k = 1}^N \epsilon(k) \left( n_{k, \downarrow} + n_{k, \uparrow} \right) + \frac{U}{N} \sum_{p, q, k = 1}^N c^\dagger_{p - k, \uparrow} c^\dagger_{q + k, \downarrow} c_{q, \downarrow} c_{p, \uparrow} $$ -where $\epsilon(k) = -2 t \cos(\frac{2 \pi k}{N})$ and $c_k = c_{k + N}$. Below is a plot of the bond dimension of the MPO produced by ITensors' default algorithm, [Renormalizer](https://github.com/shuaigroup/Renormalizer) which uses the [bipartite-graph algorithm](https://doi.org/10.1063/5.0018149), and `ITensorMPOConstruction`. +where $\epsilon(k) = -2 t \cos(\frac{2 \pi k}{N})$ and $c_k = c_{k + N}$. The code to construct the `OpSum` is shown below. -![](./docs/plot-generators/fh.png) +https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/fermi-hubbard.jl#L26-L49 -Of note is that the bond dimension of the MPO produced by Renormalizer scales as $O(N^2)$, both ITensors and ITensorMPOConstruction however produce an MPO with a bond dimension that scales as $O(N)$. Below is a table of the time it took to construct the MPO for various number of sites. Some warm up was done for the Julia calculations to avoid measuring compilation overhead. Data recorded on 2021 MacBook Pro with the M1 Max CPU and 32GB of memory. +Unlike the previous example, some more involved changes will be required to use ITensorMPOConstruction. This is because the `OpSum` has multiple operators acting on the same site, violating constraint #3. For example when $k = 0$ in the second loop we have terms of the form $c^\dagger_{p, \uparrow} c^\dagger_{q, \downarrow} c_{q, \downarrow} c_{p, \uparrow}$. You could always create a special case for $k = 0$ and rewrite it as $n_{p, \uparrow} n_{q, \downarrow}$. However if using "Electron" sites then you would also need to consider other cases such as when $p = q$, this would introduce a lot of extraneous complication. Luckily ITensorMPOConstruction provides a method to automatically perform these transformations. If you provide a set of operators for each site to `MPO_new` it will attempt to express the operators acting on each site as a single one of these "basis" operators. The code to do this is shown below. -| $N$ | ITensors | Renormalizer | ITensorMPOConstruction | -|-----|----------|--------------|------------------------| -| 10 | 0.35s | 0.26 | 0.03s | -| 20 | 27s | 3.4s | 0.15s | -| 30 | N/A | 17s | 0.41s | -| 40 | N/A | 59s | 0.93s | -| 50 | N/A | 244s | 1.8s | -| 100 | N/A | N/A | 20s | -| 200 | N/A | N/A | 310s | +https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/fermi-hubbard.jl#L51-L76 +### `OpIDSum` -The code for this example can be found in [examples/fermi-hubbard.jl](https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/fermi-hubbard.jl). Just like with ITensors, the terms in the Hamiltonian are put into an `OpSum` and then the `OpSum` is transformed into an `MPO`. The code to produce the `OpSum` is +For $N = 200$ constructing the `OpSum` takes 42s and constructing the MPO from the `OpSum` with ITensorMPOConstruction takes another 306s. For some systems constructing the `OpSum` can actually be the bottleneck. In these cases you can construct an `OpIDSum` instead. -```julia -os = OpSum{Float64}() -for k in 1:N - epsilon = -2 * t * cospi(2 * k / N) - os .+= epsilon, "Nup", k - os .+= epsilon, "Ndn", k -end - -for p in 1:N - for q in 1:N - for k in 1:N - os .+= U / N, "Cdagup", mod1(p - k, N), "Cdagdn", mod1(q + k, N), "Cdn", q, "Cup", p - end - end -end -``` +`OpIDSum` plays the same roll as `OpSum` but in a much more efficient manner. To specify an operator in a term of an `OpSum` you specify a string (the operator's name) and a site index, whereas to specify an operator in a term of an `OpIDSum` you specify an `OpID` which contains an operator index and a site. The operator index is the index of the operator in the provided basis for the site. -The astute reader will notice that this `OpSum` has multiple operators acting on the same site. For example, it contains the term `"Cdagup", 1, "Cdagdn", 1, "Cdn", 1, "Cup", 1`. If we try and pass this `OpSum` to `MPO_New` directly it will throw an error. However, `"Cdagup", 1, "Cdagdn", 1, "Cdn", 1, "Cup", 1` is equivalent to `"Nup * Ndn", 1` and by passing a set of basis operators to `MPO_New` we can automatically convert any product of operators acting on a single site into one of these single basis operators. This is accomplished by the following +For $N = 200$ constructing an `OpIDSum` takes only 0.4s. Shown below is code to construct the Hamiltonian using an `OpIDSum`. -```julia -sites = siteinds("Electron", N; conserve_qns=true) - -operatorNames = [ - "I", - "Cdn", - "Cup", - "Cdagdn", - "Cdagup", - "Ndn", - "Nup", - "Cup * Cdn", - "Cup * Cdagdn", - "Cup * Ndn", - "Cdagup * Cdn", - "Cdagup * Cdagdn", - "Cdagup * Ndn", - "Nup * Cdn", - "Nup * Cdagdn", - "Nup * Ndn", -] - -opCacheVec = [ - [OpInfo(ITensors.Op(name, n), sites[n]) for name in operatorNames] for - n in eachindex(sites) -] - -return MPO_new(os, sites; basisOpCacheVec=opCacheVec) -``` +https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/fermi-hubbard.jl#L79-L130 + +## Benchmarks: Fermi-Hubbard Hamiltonian in Momentum Space + +Below is a plot of the bond dimension of the MPO produced by ITensors' default algorithm, [Renormalizer](https://github.com/shuaigroup/Renormalizer) which uses the [bipartite-graph algorithm](https://doi.org/10.1063/5.0018149), and ITensorMPOConstruction. + +![](./docs/plot-generators/fh.png) + +Of note is that the bond dimension of the MPO produced by Renormalizer scales as $O(N^2)$, both ITensors and ITensorMPOConstruction however produce an MPO with a bond dimension that scales as $O(N)$. + +Below is a table of the time it took to construct the MPO (including the time it took to specify the Hamiltonian) for various number of sites. For ITensorMPOConstruction an `OpIDSum` was used. Some warm up was done for the Julia calculations to avoid measuring compilation overhead. + +| $N$ | ITensors | Renormalizer | ITensorMPOConstruction | +|-----|----------|--------------|------------------------| +| 10 | 0.35s | 0.26 | 0.04s | +| 20 | 27s | 3.4s | 0.10s | +| 30 | N/A | 17s | 0.24s | +| 40 | N/A | 59s | 0.59s | +| 50 | N/A | 244s | 1.2s | +| 100 | N/A | N/A | 16s | +| 200 | N/A | N/A | 283s | -## Examples: Electronic Structure Hamiltonian +## Benchmarks: Electronic Structure Hamiltonian -After looking at the previous example one might assume that the relative speed of ITensorMPOConstruction over Renormalizer might be due to the fact that for the Fermi-Hubbard Hamiltonian ITensorMPOConstruction is able to construct a much more compact MPO. In the case of the electronic structure Hamiltonian all algorithms produce MPOs of similar bond dimensions. +After looking at the previous example you might assume that the relative speed of ITensorMPOConstruction over Renormalizer might be due to the fact that for the Fermi-Hubbard Hamiltonian ITensorMPOConstruction is able to construct a more compact MPO. In the case of the electronic structure Hamiltonian all algorithms produce MPOs of similar bond dimensions. ![](./docs/plot-generators/es.png) -However ITensorMPOConstruction is still an order of magnitude faster than Renormalizer. The code for this example can be found in [examples/electronic-structure.jl](https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/electronic-structure.jl). The run time to generate these MPOs are shown below (again on a 2021 MacBook Pro with the M1 Max CPU and 32GB of memory). +However, ITensorMPOConstruction is still an order of magnitude faster than Renormalizer. The code for this example can be found in [examples/electronic-structure.jl](https://github.com/ITensor/ITensorMPOConstruction.jl/blob/main/examples/electronic-structure.jl). The run time to generate these MPOs (including the time it took to specify the Hamiltonian) are shown below. For ITensorMPOConstruction an `OpIDSum` was used. | $N$ | ITensors | Renormalizer | ITensorMPOConstruction | |-----|----------|--------------|------------------------| -| 10 | 3.0s | 2.1s | 0.37s | -| 20 | 498s | 61s | 7.1s | -| 30 | N/A | 458s | 46s | -| 40 | N/A | 2250s | 183 | -| 50 | N/A | N/A | 510s | +| 10 | 3.0s | 2.1s | 0.31s | +| 20 | 498s | 61s | 5.9s | +| 30 | N/A | 458s | 36s | +| 40 | N/A | 2250s | 162s | +| 50 | N/A | N/A | 504s | +| 60 | N/A | N/A | 1323s | diff --git a/examples/electronic-structure.jl b/examples/electronic-structure.jl index cbc3018..435d68c 100644 --- a/examples/electronic-structure.jl +++ b/examples/electronic-structure.jl @@ -1,41 +1,43 @@ using ITensorMPOConstruction using ITensors -function electronic_structure_example( +function electronic_structure( N::Int, complexBasisFunctions::Bool; useITensorsAlg::Bool=false )::MPO - coeff() = rand() + 1im * complexBasisFunctions * rand() + coeff() = !complexBasisFunctions ? rand() : rand() + 1im * rand() os = complexBasisFunctions ? OpSum{ComplexF64}() : OpSum{Float64}() - for a in 1:N - for b in a:N - epsilon_ab = coeff() - os .+= epsilon_ab, "Cdagup", a, "Cup", b - os .+= epsilon_ab, "Cdagdn", a, "Cdn", b - - a == b && continue - os .+= conj(epsilon_ab), "Cdagup", b, "Cup", a - os .+= conj(epsilon_ab), "Cdagdn", b, "Cdn", a + @time "\tConstructing OpSum" let + for a in 1:N + for b in a:N + epsilon_ab = coeff() + os .+= epsilon_ab, "Cdagup", a, "Cup", b + os .+= epsilon_ab, "Cdagdn", a, "Cdn", b + + a == b && continue + os .+= conj(epsilon_ab), "Cdagup", b, "Cup", a + os .+= conj(epsilon_ab), "Cdagdn", b, "Cdn", a + end end - end - for j in 1:N - for s_j in ("dn", "up") - for k in 1:N - s_k = s_j - (s_k, k) <= (s_j, j) && continue + for j in 1:N + for s_j in ("dn", "up") + for k in 1:N + s_k = s_j + (s_k, k) <= (s_j, j) && continue - for l in 1:N - for s_l in ("dn", "up") - (s_l, l) <= (s_j, j) && continue + for l in 1:N + for s_l in ("dn", "up") + (s_l, l) <= (s_j, j) && continue - for m in 1:N - s_m = s_l - (s_m, m) <= (s_k, k) && continue + for m in 1:N + s_m = s_l + (s_m, m) <= (s_k, k) && continue - value = coeff() - os .+= value, "Cdag$s_j", j, "Cdag$s_l", l, "C$s_m", m, "C$s_k", k - os .+= conj(value), "Cdag$s_k", k, "Cdag$s_m", m, "C$s_l", l, "C$s_j", j + value = coeff() + os .+= value, "Cdag$s_j", j, "Cdag$s_l", l, "C$s_m", m, "C$s_k", k + os .+= conj(value), "Cdag$s_k", k, "Cdag$s_m", m, "C$s_l", l, "C$s_j", j + end end end end @@ -47,7 +49,7 @@ function electronic_structure_example( ## The only additional step required is to provide an operator basis in which to represent the OpSum. if useITensorsAlg - return MPO(os, sites) + return @time "\tConstrucing MPO" MPO(os, sites) else operatorNames = [ "I", @@ -73,18 +75,102 @@ function electronic_structure_example( n in eachindex(sites) ] - return MPO_new(os, sites; basisOpCacheVec=opCacheVec) + return @time "\tConstrucing MPO" MPO_new(os, sites; basisOpCacheVec=opCacheVec) end end -let - N = 25 +function electronic_structure_OpIDSum(N::Int, complexBasisFunctions::Bool)::MPO + operatorNames = [ + [ + "I", + "Cdn", + "Cup", + "Cdagdn", + "Cdagup", + "Ndn", + "Nup", + "Cup * Cdn", + "Cup * Cdagdn", + "Cup * Ndn", + "Cdagup * Cdn", + "Cdagup * Cdagdn", + "Cdagup * Ndn", + "Nup * Cdn", + "Nup * Cdagdn", + "Nup * Ndn", + ] for _ in 1:N + ] + + ↓ = false + ↑ = true + + opC(k::Int, spin::Bool) = OpID(2 + spin, k) + opCdag(k::Int, spin::Bool) = OpID(4 + spin, k) + opN(k::Int, spin::Bool) = OpID(6 + spin, k) + + coeff() = !complexBasisFunctions ? rand() : rand() + 1im * rand() - # Ensure compilation - electronic_structure_example(5, false) + os = complexBasisFunctions ? OpIDSum{ComplexF64}() : OpIDSum{Float64}() + @time "\tConstructing OpIDSum" let + for a in 1:N + for b in a:N + epsilon_ab = coeff() + push!(os, epsilon_ab, opCdag(a, ↑), opC(b, ↑)) + push!(os, epsilon_ab, opCdag(a, ↓), opC(b, ↓)) + + a == b && continue + push!(os, conj(epsilon_ab), opCdag(b, ↑), opC(a, ↑)) + push!(os, conj(epsilon_ab), opCdag(b, ↓), opC(a, ↓)) + end + end + + for j in 1:N + for s_j in (↓, ↑) + for k in 1:N + s_k = s_j + (s_k, k) <= (s_j, j) && continue + + for l in 1:N + for s_l in (↓, ↑) + (s_l, l) <= (s_j, j) && continue + + for m in 1:N + s_m = s_l + (s_m, m) <= (s_k, k) && continue + + value = coeff() + push!(os, value, opCdag(j, s_j), opCdag(l, s_l), opC(m, s_m), opC(k, s_k)) + push!( + os, conj(value), opCdag(k, s_k), opCdag(m, s_m), opC(l, s_l), opC(j, s_j) + ) + end + end + end + end + end + end + end + + sites = siteinds("Electron", N; conserve_qns=true) + return @time "\tConstructing MPO" MPO_new( + os, sites, operatorNames; basisOpCacheVec=operatorNames + ) +end + +let + N = 20 + useITensorsAlg = false println("Constructing the Electronic Structure MPO for $N orbitals") - @time mpo = electronic_structure_example(N, false) + @time "Total" mpo = electronic_structure(N, false; useITensorsAlg=useITensorsAlg) + println("The maximum bond dimension is $(maxlinkdim(mpo))") +end + +let + N = 20 + + println("Constructing the Electronic Structure MPO for $N orbitals using OpIDSum") + @time "Total" mpo = electronic_structure_OpIDSum(N, false) println("The maximum bond dimension is $(maxlinkdim(mpo))") end diff --git a/examples/fermi-hubbard.jl b/examples/fermi-hubbard.jl index d56b88d..9f12940 100644 --- a/examples/fermi-hubbard.jl +++ b/examples/fermi-hubbard.jl @@ -1,32 +1,84 @@ using ITensorMPOConstruction using ITensors -function Fermi_Hubbard_momentum_space( +function Fermi_Hubbard_real_space( N::Int, t::Real=1.0, U::Real=4.0; useITensorsAlg::Bool=false )::MPO - ## Create the OpSum as per usual. os = OpSum{Float64}() - for k in 1:N - epsilon = -2 * t * cospi(2 * k / N) - os .+= epsilon, "Nup", k - os .+= epsilon, "Ndn", k + for i in 1:N + for j in [mod1(i + 1, N), mod1(i - 1, N)] + os .+= -t, "Cdagup", i, "Cup", j + os .+= -t, "Cdagdn", i, "Cdn", j + end + + os .+= U, "Nup * Ndn", i + end + + sites = siteinds("Electron", N; conserve_qns=true) + + if useITensorsAlg + return MPO(os, sites) + else + return MPO_new(os, sites) end +end - for p in 1:N - for q in 1:N - for k in 1:N - os .+= U / N, "Cdagup", mod1(p - k, N), "Cdagdn", mod1(q + k, N), "Cdn", q, "Cup", p +function Fermi_Hubbard_momentum_space( + N::Int, t::Real=1.0, U::Real=4.0; useITensorsAlg::Bool=false +)::MPO + os = OpSum{Float64}() + + @time "Constructing OpSum\t" let + for k in 1:N + epsilon = -2 * t * cospi(2 * k / N) + os .+= epsilon, "Nup", k + os .+= epsilon, "Ndn", k + end + + for p in 1:N + for q in 1:N + for k in 1:N + os .+= U / N, + "Cdagup", mod1(p - k, N), "Cdagdn", mod1(q + k, N), "Cdn", q, "Cup", + p + end end end end sites = siteinds("Electron", N; conserve_qns=true) - ## The only additional step required is to provide an operator basis in which to represent the OpSum. if useITensorsAlg - return MPO(os, sites) + return @time "\tConstructing MPO" MPO(os, sites) else operatorNames = [ + [ + "I", + "Cdn", + "Cup", + "Cdagdn", + "Cdagup", + "Ndn", + "Nup", + "Cup * Cdn", + "Cup * Cdagdn", + "Cup * Ndn", + "Cdagup * Cdn", + "Cdagup * Cdagdn", + "Cdagup * Ndn", + "Nup * Cdn", + "Nup * Cdagdn", + "Nup * Ndn", + ] for _ in 1:N + ] + + return @time "\tConstructing MPO" MPO_new(os, sites; basisOpCacheVec=operatorNames) + end +end + +function Fermi_Hubbard_momentum_space_OpIDSum(N::Int, t::Real=1.0, U::Real=4.0)::MPO + operatorNames = [ + [ "I", "Cdn", "Cup", @@ -43,25 +95,63 @@ function Fermi_Hubbard_momentum_space( "Nup * Cdn", "Nup * Cdagdn", "Nup * Ndn", - ] + ] for _ in 1:N + ] - opCacheVec = [ - [OpInfo(ITensors.Op(name, n), sites[n]) for name in operatorNames] for - n in eachindex(sites) - ] + ↓ = false + ↑ = true - return MPO_new(os, sites; basisOpCacheVec=opCacheVec) + opC(k::Int, spin::Bool) = OpID(2 + spin, mod1(k, N)) + opCdag(k::Int, spin::Bool) = OpID(4 + spin, mod1(k, N)) + opN(k::Int, spin::Bool) = OpID(6 + spin, mod1(k, N)) + + os = OpIDSum{Float64}() + + @time "\tConstructing OpIDSum" let + for k in 1:N + epsilon = -2 * t * cospi(2 * k / N) + push!(os, epsilon, opN(k, ↑)) + push!(os, epsilon, opN(k, ↓)) + end + + for p in 1:N + for q in 1:N + for k in 1:N + push!(os, U / N, opCdag(p - k, ↑), opCdag(q + k, ↓), opC(q, ↓), opC(p, ↑)) + end + end + end end + + sites = siteinds("Electron", N; conserve_qns=true) + return @time "\tConstructing MPO" MPO_new( + os, sites, operatorNames; basisOpCacheVec=operatorNames + ) end let N = 50 + useITensorsAlg = false - # Ensure compilation - Fermi_Hubbard_momentum_space(10) + println("Constructing the Fermi-Hubbard real space MPO for $N sites") + mpo = Fermi_Hubbard_real_space(N; useITensorsAlg=useITensorsAlg) + println("The maximum bond dimension is $(maxlinkdim(mpo))") +end + +let + N = 50 + useITensorsAlg = false println("Constructing the Fermi-Hubbard momentum space MPO for $N sites") - @time mpo = Fermi_Hubbard_momentum_space(N) + mpo = Fermi_Hubbard_momentum_space(N; useITensorsAlg=useITensorsAlg) + println("The maximum bond dimension is $(maxlinkdim(mpo))") +end + +let + N = 50 + + println("Constructing the Fermi-Hubbard momentum space MPO for $N sites using OpIDSum") + mpo = Fermi_Hubbard_momentum_space_OpIDSum(N) println("The maximum bond dimension is $(maxlinkdim(mpo))") end diff --git a/examples/haldane-shastry.jl b/examples/haldane-shastry.jl new file mode 100644 index 0000000..abd7b97 --- /dev/null +++ b/examples/haldane-shastry.jl @@ -0,0 +1,37 @@ +using ITensorMPOConstruction +using ITensors + +function halden_shastry_mpo_from_OpSum(N::Int, J::Real; useITensorsAlg::Bool=false)::MPO + os = OpSum() + @time "Constructing OpSum" for n in 1:N + for m in (n + 1):N + factor = J * π^2 / N^2 / sinpi((n - m) / N)^2 + os .+= factor, "Sz", n, "Sz", m + os .+= factor / 2, "S+", n, "S-", m + os .+= factor / 2, "S-", n, "S+", m + end + end + + sites = siteinds("S=1/2", N; conserve_sz=true) + + if useITensorsAlg + return @time "Constructing MPO" MPO(os, sites) + else + return @time "Constructing MPO" MPO_new(os, sites) + end +end + +# function ground_state_energy(N::Int, Mtimes2::Int, J::Real) +# M = abs(Mtimes2) / 2 +# return J * π^2 * ((N - 1 / N) / 24 + M * ((M^2 - 1) / (3 * N^2) - 1 / 4)) +# end + +let + N = 800 + J = 1.0 + + H, sites = halden_shastry_mpo_from_OpSum(N, 1.0) + @show maxlinkdim(H) +end + +nothing; diff --git a/examples/intro.jl b/examples/intro.jl new file mode 100644 index 0000000..1a76326 --- /dev/null +++ b/examples/intro.jl @@ -0,0 +1,36 @@ +using ITensorMPOConstruction +using ITensors + +function foo(N::Int; useITensorsAlg::Bool=false)::MPO + os = OpSum{Float64}() + @time "Constructing OpSum" for i in 1:N + for j in (i + 1):N + for k in (j + 1):N + for l in (k + 1):N + os .+= 1, "X", i, "X", j, "X", k, "X", l + os .+= 1, "Y", i, "Y", j, "Y", k, "Y", l + os .+= 1, "Z", i, "Z", j, "Z", k, "Z", l + end + end + end + end + + sites = siteinds("Qubit", N) + + if useITensorsAlg + return @time "Constructing MPO" MPO(os, sites) + else + return @time "Constructing MPO" MPO_new(os, sites) + end +end + +let + for N in 5:5:100 + println("N = $N") + H = foo(N; useITensorsAlg=false) + @show maxlinkdim(H) + println() + end +end + +nothing; diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl index 7495b99..7954f24 100644 --- a/src/MPOConstruction.jl +++ b/src/MPOConstruction.jl @@ -53,16 +53,17 @@ function for_non_zeros_batch(f::Function, A::SparseMatrixCSC, cols::Vector{Int}) end @timeit function at_site!( + ValType::Type{<:Number}, graphs::Dict{QN,MPOGraph{C}}, n::Int, sites::Vector{<:Index}, tol::Real, opCacheVec::OpCacheVec, -)::Tuple{Dict{QN,MPOGraph{C}},BlockSparseMatrix{C},Index} where {C} +)::Tuple{Dict{QN,MPOGraph{C}},BlockSparseMatrix{ValType},Index} where {C} hasQNs = hasqns(sites) nextGraphs = Dict{QN,MPOGraph{C}}() - matrix = BlockSparseMatrix{C}() + matrix = BlockSparseMatrix{ValType}() qi = Vector{Pair{QN,Int}}() outgoingLinkOffset = 0 @@ -194,7 +195,9 @@ end llinks = Vector{Index}(undef, N + 1) for n in 0:N - graphs, symbolicMatrix, llinks[n + 1] = at_site!(graphs, n, sites, tol, opCacheVec) + graphs, symbolicMatrix, llinks[n + 1] = at_site!( + ValType, graphs, n, sites, tol, opCacheVec + ) # For the 0th iteration we only care about constructing the graphs for the next site. n == 0 && continue @@ -239,27 +242,46 @@ end return H end -function svdMPO_new(os::OpIDSum{C}, opCacheVec::OpCacheVec, sites; kwargs...)::MPO where {C} - # Function barrier to improve type stability - ValType = determine_val_type(os) - return svdMPO_new(ValType, os, opCacheVec, sites; kwargs...) -end - function MPO_new( + ValType::Type{<:Number}, os::OpIDSum, sites::Vector{<:Index}, opCacheVec::OpCacheVec; - basisOpCacheVec::Union{Nothing,OpCacheVec}=nothing, - kwargs..., + basisOpCacheVec=nothing, + tol::Real=-1, )::MPO + opCacheVec = to_OpCacheVec(sites, opCacheVec) + basisOpCacheVec = to_OpCacheVec(sites, basisOpCacheVec) os, opCacheVec = prepare_opID_sum!(os, sites, opCacheVec, basisOpCacheVec) + return svdMPO_new(ValType, os, opCacheVec, sites; tol=tol) +end + +function MPO_new( + os::OpIDSum, sites::Vector{<:Index}, opCacheVec; basisOpCacheVec=nothing, tol::Real=-1 +)::MPO + opCacheVec = to_OpCacheVec(sites, opCacheVec) + ValType = determine_val_type(os, opCacheVec) + return MPO_new(ValType, os, sites, opCacheVec; basisOpCacheVec=basisOpCacheVec, tol=tol) +end - return svdMPO_new(os, opCacheVec, sites; kwargs...) +function MPO_new( + ValType::Type{<:Number}, + os::OpSum, + sites::Vector{<:Index}; + basisOpCacheVec=nothing, + tol::Real=-1, +)::MPO + opIDSum, opCacheVec = op_sum_to_opID_sum(os, sites) + return MPO_new( + ValType, opIDSum, sites, opCacheVec; basisOpCacheVec=basisOpCacheVec, tol=tol + ) end -function MPO_new(os::OpSum, sites::Vector{<:Index}; kwargs...)::MPO +function MPO_new( + os::OpSum, sites::Vector{<:Index}; basisOpCacheVec=nothing, tol::Real=-1 +)::MPO opIDSum, opCacheVec = op_sum_to_opID_sum(os, sites) - return MPO_new(opIDSum, sites, opCacheVec; kwargs...) + return MPO_new(opIDSum, sites, opCacheVec; basisOpCacheVec=basisOpCacheVec, tol=tol) end function sparsity(mpo::MPO)::Float64 diff --git a/src/OpIDSum.jl b/src/OpIDSum.jl index 8a3202c..ba527dc 100644 --- a/src/OpIDSum.jl +++ b/src/OpIDSum.jl @@ -14,6 +14,26 @@ end OpCacheVec = Vector{Vector{OpInfo}} +function to_OpCacheVec(sites::Vector{<:Index}, ops::OpCacheVec)::OpCacheVec + length(sites) != length(ops) && + error("Mismatch in the number of sites in `sites` and `ops`.") + any(ops[i][1].matrix != I for i in 1:length(sites)) && + error("The first operator on each site must be the identity.") + return ops +end + +function to_OpCacheVec(sites::Vector{<:Index}, ops::Vector{Vector{String}})::OpCacheVec + length(sites) != length(ops) && + error("Mismatch in the number of sites in `sites` and `ops`.") + any(ops[i][1] != "I" for i in 1:length(sites)) && + error("The first operator on each site must be the identity.") + return [[OpInfo(ITensors.Op(op, n), sites[n]) for op in ops[n]] for n in 1:length(sites)] +end + +function to_OpCacheVec(sites::Vector{<:Index}, ::Nothing)::Nothing + return nothing +end + struct OpID id::Int16 n::Int16 @@ -71,13 +91,10 @@ function add_to_scalar!(os::OpIDSum{C}, i::Integer, scalar::C)::Nothing where {C return nothing end -# TODO: Define as `C`. Rename `coefficient_type`. -function determine_val_type(os::OpIDSum{C}) where {C} - for i in eachindex(os) - scalar, ops = os[i] - (!isreal(scalar)) && return ComplexF64 - end - +function determine_val_type(os::OpIDSum{C}, opCacheVec::OpCacheVec) where {C} + !all(isreal(scalar) for scalar in os.scalars) && return ComplexF64 + !all(isreal(op.matrix) for opsOfSite in opCacheVec for op in opsOfSite) && + return ComplexF64 return Float64 end