Skip to content

Commit 3cc10c1

Browse files
mtfishmanmtfishman
authored
Add solver examples and export contract (#20)
* Add solver examples and export contract * Bump to v0.2.3
1 parent 39f1a9b commit 3cc10c1

14 files changed

+414
-41
lines changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "ITensorMPS"
22
uuid = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
33
authors = ["Matthew Fishman <[email protected]>", "Miles Stoudenmire <[email protected]>"]
4-
version = "0.2.2"
4+
version = "0.2.3"
55

66
[deps]
77
ITensorTDVP = "25707e16-a4db-4a07-99d9-4d67b7af0342"

‎examples/01_tdvp.jl‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
using ITensorMPS: MPO, OpSum, dmrg, inner, random_mps, siteinds, tdvp
2+
3+
function main()
4+
n = 10
5+
s = siteinds("S=1/2", n)
6+
7+
function heisenberg(n)
8+
os = OpSum()
9+
for j in 1:(n - 1)
10+
os += 0.5, "S+", j, "S-", j + 1
11+
os += 0.5, "S-", j, "S+", j + 1
12+
os += "Sz", j, "Sz", j + 1
13+
end
14+
return os
15+
end
16+
17+
H = MPO(heisenberg(n), s)
18+
ψ = random_mps(s, "↑"; linkdims=10)
19+
20+
@show inner(ψ', H, ψ) / inner(ψ, ψ)
21+
22+
ϕ = tdvp(
23+
H,
24+
-20.0,
25+
ψ;
26+
time_step=-1.0,
27+
maxdim=30,
28+
cutoff=1e-10,
29+
normalize=true,
30+
reverse_step=false,
31+
outputlevel=1,
32+
)
33+
@show inner(ϕ', H, ϕ) / inner(ϕ, ϕ)
34+
35+
e2, ϕ2 = dmrg(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10)
36+
@show inner(ϕ2', H, ϕ2) / inner(ϕ2, ϕ2), e2
37+
38+
return nothing
39+
end
40+
41+
main()

‎examples/02_dmrg-x.jl‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
using ITensorMPS: MPO, MPS, OpSum, dmrg_x, inner, siteinds
2+
using Random: Random
3+
4+
function main()
5+
function heisenberg(n; h=zeros(n))
6+
os = OpSum()
7+
for j in 1:(n - 1)
8+
os += 0.5, "S+", j, "S-", j + 1
9+
os += 0.5, "S-", j, "S+", j + 1
10+
os += "Sz", j, "Sz", j + 1
11+
end
12+
for j in 1:n
13+
if h[j] ≠ 0
14+
os -= h[j], "Sz", j
15+
end
16+
end
17+
return os
18+
end
19+
20+
n = 10
21+
s = siteinds("S=1/2", n)
22+
23+
Random.seed!(12)
24+
25+
# MBL when W > 3.5-4
26+
W = 12
27+
# Random fields h ∈ [-W, W]
28+
h = W * (2 * rand(n) .- 1)
29+
H = MPO(heisenberg(n; h), s)
30+
31+
initstate = rand(["↑", "↓"], n)
32+
ψ = MPS(s, initstate)
33+
e, ϕ = dmrg_x(H, ψ; nsweeps=10, maxdim=20, cutoff=1e-10, normalize=true, outputlevel=1)
34+
35+
@show inner(ψ', H, ψ) / inner(ψ, ψ)
36+
@show inner(H, ψ, H, ψ) - inner(ψ', H, ψ)^2
37+
@show inner(ϕ', H, ϕ) / inner(ϕ, ϕ), e
38+
@show inner(H, ϕ, H, ϕ) - inner(ϕ', H, ϕ)^2
39+
return nothing
40+
end
41+
42+
main()

‎examples/03_models.jl‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
using ITensorMPS: OpSum
2+
3+
function heisenberg(n; J=1.0, J2=0.0)
4+
ℋ = OpSum()
5+
if !iszero(J)
6+
for j in 1:(n - 1)
7+
ℋ += J / 2, "S+", j, "S-", j + 1
8+
ℋ += J / 2, "S-", j, "S+", j + 1
9+
ℋ += J, "Sz", j, "Sz", j + 1
10+
end
11+
end
12+
if !iszero(J2)
13+
for j in 1:(n - 2)
14+
ℋ += J2 / 2, "S+", j, "S-", j + 2
15+
ℋ += J2 / 2, "S-", j, "S+", j + 2
16+
ℋ += J2, "Sz", j, "Sz", j + 2
17+
end
18+
end
19+
return ℋ
20+
end

‎examples/03_tdvp_time_dependent.jl‎

Lines changed: 153 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,153 @@
1+
using ITensors: @disable_warn_order
2+
using ITensorMPS: MPO, MPS, contract, inner, random_mps, siteinds, tdvp
3+
using LinearAlgebra: norm
4+
using Random: Random
5+
6+
include("03_models.jl")
7+
include("03_updaters.jl")
8+
9+
function main()
10+
Random.seed!(1234)
11+
12+
# Time dependent Hamiltonian is:
13+
# H(t) = H₁(t) + H₂(t) + …
14+
# = f₁(t) H₁(0) + f₂(t) H₂(0) + …
15+
# = cos(ω₁t) H₁(0) + cos(ω₂t) H₂(0) + …
16+
17+
# Number of sites
18+
n = 6
19+
20+
# How much information to output from TDVP
21+
# Set to 2 to get information about each bond/site
22+
# evolution, and 3 to get information about the
23+
# updater.
24+
outputlevel = 3
25+
26+
# Frequency of time dependent terms
27+
ω₁ = 0.1
28+
ω₂ = 0.2
29+
30+
# Nearest and next-nearest neighbor
31+
# Heisenberg couplings.
32+
J₁ = 1.0
33+
J₂ = 1.0
34+
35+
time_step = 0.1
36+
time_stop = 1.0
37+
38+
# nsite-update TDVP
39+
nsite = 2
40+
41+
# Starting state bond/link dimension.
42+
# A product state starting state can
43+
# cause issues for TDVP without
44+
# subspace expansion.
45+
start_linkdim = 4
46+
47+
# TDVP truncation parameters
48+
maxdim = 100
49+
cutoff = 1e-8
50+
51+
tol = 1e-15
52+
53+
@show n
54+
@show ω₁, ω₂
55+
@show J₁, J₂
56+
@show maxdim, cutoff, nsite
57+
@show start_linkdim
58+
@show time_step, time_stop
59+
60+
ω⃗ = (ω₁, ω₂)
61+
f⃗ = map(ω -> (t -> cos(ω * t)), ω⃗)
62+
63+
# H₀ = H(0) = H₁(0) + H₂(0) + …
64+
ℋ₁₀ = heisenberg(n; J=J₁, J2=0.0)
65+
ℋ₂₀ = heisenberg(n; J=0.0, J2=J₂)
66+
ℋ⃗₀ = (ℋ₁₀, ℋ₂₀)
67+
68+
s = siteinds("S=1/2", n)
69+
70+
H⃗₀ = map(ℋ₀ -> MPO(ℋ₀, s), ℋ⃗₀)
71+
72+
# Initial state, ψ₀ = ψ(0)
73+
# Initialize as complex since that is what OrdinaryDiffEq.jl/DifferentialEquations.jl
74+
# expects.
75+
ψ₀ = complex.(random_mps(s, j -> isodd(j) ? "↑" : "↓"; linkdims=start_linkdim))
76+
77+
@show norm(ψ₀)
78+
79+
println()
80+
println("#"^100)
81+
println("Running TDVP with ODE updater")
82+
println("#"^100)
83+
println()
84+
85+
ψₜ_ode = tdvp(
86+
-im * TimeDependentSum(f⃗, H⃗₀),
87+
time_stop,
88+
ψ₀;
89+
updater=ode_updater,
90+
updater_kwargs=(; reltol=tol, abstol=tol),
91+
time_step,
92+
maxdim,
93+
cutoff,
94+
nsite,
95+
outputlevel,
96+
)
97+
98+
println()
99+
println("Finished running TDVP with ODE updater")
100+
println()
101+
102+
println()
103+
println("#"^100)
104+
println("Running TDVP with Krylov updater")
105+
println("#"^100)
106+
println()
107+
108+
ψₜ_krylov = tdvp(
109+
-im * TimeDependentSum(f⃗, H⃗₀),
110+
time_stop,
111+
ψ₀;
112+
updater=krylov_updater,
113+
updater_kwargs=(; tol, eager=true),
114+
time_step,
115+
cutoff,
116+
nsite,
117+
outputlevel,
118+
)
119+
120+
println()
121+
println("Finished running TDVP with Krylov updater")
122+
println()
123+
124+
println()
125+
println("#"^100)
126+
println("Running full state evolution with ODE updater")
127+
println("#"^100)
128+
println()
129+
130+
@disable_warn_order begin
131+
ψₜ_full, _ = ode_updater(
132+
-im * TimeDependentSum(f⃗, contract.(H⃗₀)),
133+
contract(ψ₀);
134+
internal_kwargs=(; time_step=time_stop, outputlevel),
135+
reltol=tol,
136+
abstol=tol,
137+
)
138+
end
139+
140+
println()
141+
println("Finished full state evolution with ODE updater")
142+
println()
143+
144+
@show norm(ψₜ_ode)
145+
@show norm(ψₜ_krylov)
146+
@show norm(ψₜ_full)
147+
148+
@show 1 - abs(inner(contract(ψₜ_ode), ψₜ_full))
149+
@show 1 - abs(inner(contract(ψₜ_krylov), ψₜ_full))
150+
return nothing
151+
end
152+
153+
main()

‎examples/03_updaters.jl‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
using Compat: @compat
2+
using ITensors: ITensor, array, inds, itensor
3+
using ITensorMPS: TimeDependentSum, to_vec
4+
using KrylovKit: exponentiate
5+
using OrdinaryDiffEq: ODEProblem, Tsit5, solve
6+
7+
function ode_updater(operator, init; internal_kwargs, alg=Tsit5(), kwargs...)
8+
@compat (; current_time, time_step) = (; current_time=zero(Bool), internal_kwargs...)
9+
time_span = typeof(time_step).((current_time, current_time + time_step))
10+
init_vec, to_itensor = to_vec(init)
11+
f(init::ITensor, p, t) = operator(t)(init)
12+
f(init_vec::Vector, p, t) = to_vec(f(to_itensor(init_vec), p, t))[1]
13+
prob = ODEProblem(f, init_vec, time_span)
14+
sol = solve(prob, alg; kwargs...)
15+
state_vec = sol.u[end]
16+
return to_itensor(state_vec), (;)
17+
end
18+
19+
function krylov_updater(operator, init; internal_kwargs, kwargs...)
20+
@compat (; current_time, time_step) = (; current_time=zero(Bool), internal_kwargs...)
21+
state, info = exponentiate(operator(current_time), time_step, init; kwargs...)
22+
return state, (; info)
23+
end

‎examples/04_tdvp_observers.jl‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
using ITensorMPS: MPO, MPS, OpSum, expect, inner, siteinds, tdvp
2+
using Observers: observer
3+
4+
function main()
5+
function heisenberg(N)
6+
os = OpSum()
7+
for j in 1:(N - 1)
8+
os += 0.5, "S+", j, "S-", j + 1
9+
os += 0.5, "S-", j, "S+", j + 1
10+
os += "Sz", j, "Sz", j + 1
11+
end
12+
return os
13+
end
14+
15+
N = 10
16+
s = siteinds("S=1/2", N; conserve_qns=true)
17+
H = MPO(heisenberg(N), s)
18+
19+
step(; sweep) = sweep
20+
current_time(; current_time) = current_time
21+
return_state(; state) = state
22+
measure_sz(; state) = expect(state, "Sz"; sites=length(state) ÷ 2)
23+
obs = observer(
24+
"steps" => step, "times" => current_time, "states" => return_state, "sz" => measure_sz
25+
)
26+
27+
init = MPS(s, n -> isodd(n) ? "Up" : "Dn")
28+
state = tdvp(
29+
H, -1.0im, init; time_step=-0.1im, cutoff=1e-12, (step_observer!)=obs, outputlevel=1
30+
)
31+
32+
println("\nResults")
33+
println("=======")
34+
for n in 1:length(obs.steps)
35+
print("step = ", obs.steps[n])
36+
print(", time = ", round(obs.times[n]; digits=3))
37+
print(", |⟨ψⁿ|ψⁱ⟩| = ", round(abs(inner(obs.states[n], init)); digits=3))
38+
print(", |⟨ψⁿ|ψᶠ⟩| = ", round(abs(inner(obs.states[n], state)); digits=3))
39+
print(", ⟨Sᶻ⟩ = ", round(obs.sz[n]; digits=3))
40+
println()
41+
end
42+
return nothing
43+
end
44+
45+
main()

‎examples/Project.toml‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
[deps]
2+
Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
3+
ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
4+
ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
5+
KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
6+
Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0"
7+
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
8+
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
9+
10+
[compat]
11+
ITensors = "0.6.7"

‎src/ITensorMPS.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ using .Experimental: Experimental
1010
include("Deprecated.jl")
1111
using .Deprecated: Deprecated, dmrg
1212
export dmrg
13+
@reexport using ITensors: contract
1314
@reexport using ITensors.ITensorMPS:
1415
@OpName_str,
1516
@SiteType_str,

‎test/Project.toml‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,12 @@
11
[deps]
2+
Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
23
ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
34
ITensorTDVP = "25707e16-a4db-4a07-99d9-4d67b7af0342"
45
ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
6+
KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
7+
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
8+
Observers = "338f10d5-c7f1-4033-a7d1-f9dec39bcaa0"
9+
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
10+
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
11+
Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
512
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

0 commit comments

Comments
 (0)