|
1 | 1 | @eval module $(gensym()) |
2 | | -using DataGraphs: vertex_data |
3 | | -using Dictionaries: Dictionary, getindices |
4 | | -using Graphs: add_edge!, add_vertex!, rem_edge!, rem_vertex!, vertices |
5 | | -using ITensorMPS: ITensorMPS |
6 | | -using ITensorNetworks.ITensorsExtensions: replace_vertices |
7 | | -using ITensorNetworks.ModelHamiltonians: ModelHamiltonians |
| 2 | +using Graphs: vertices |
8 | 3 | using ITensorNetworks: ITensorNetworks, OpSum, siteinds, ttn |
9 | | -using ITensors.NDTensors: matrix |
10 | | -using ITensors: ITensors, @disable_warn_order, ITensor, Index, combinedind, combiner, |
11 | | - contract, dag, inds, removeqns |
12 | | -using KrylovKit: eigsolve |
13 | | -using LinearAlgebra: eigvals, norm |
14 | | -using NamedGraphs.GraphsExtensions: leaf_vertices, post_order_dfs_vertices |
15 | | -using NamedGraphs.NamedGraphGenerators: named_comb_tree, named_grid |
16 | | -using Test: @test, @test_broken, @testset |
17 | | - |
18 | | -function to_matrix(t::ITensor) |
19 | | - c = combiner(inds(t; plev = 0)) |
20 | | - tc = (t * c) * dag(c') |
21 | | - cind = combinedind(c) |
22 | | - return matrix(tc, cind', cind) |
23 | | -end |
| 4 | +using ITensors: ITensors |
| 5 | +using NamedGraphs.NamedGraphGenerators: named_grid |
| 6 | +using Test: @test, @testset |
24 | 7 |
|
25 | 8 | @testset "OpSum to TTN converter" begin |
26 | | - @testset "OpSum to TTN" begin |
27 | | - # small comb tree |
28 | | - auto_fermion_enabled = ITensors.using_auto_fermion() |
29 | | - tooth_lengths = fill(2, 3) |
30 | | - c = named_comb_tree(tooth_lengths) |
31 | | - |
32 | | - is = siteinds("S=1/2", c) |
33 | | - |
34 | | - # linearized version |
35 | | - linear_order = [4, 1, 2, 5, 3, 6] |
36 | | - vmap = Dictionary(collect(vertices(is))[linear_order], eachindex(linear_order)) |
37 | | - sites = only.(collect(vertex_data(is)))[linear_order] |
38 | | - |
39 | | - # test with next-to-nearest-neighbor Ising Hamiltonian |
40 | | - J1 = -1 |
41 | | - J2 = 2 |
42 | | - h = 0.5 |
43 | | - H = ModelHamiltonians.ising(c; J1 = J1, J2 = J2, h = h) |
44 | | - # add combination of longer range interactions |
45 | | - Hlr = copy(H) |
46 | | - Hlr += 5, "Z", (1, 2), "Z", (2, 2), "Z", (3, 2) |
47 | | - Hlr += -4, "Z", (1, 1), "Z", (2, 2), "Z", (3, 1) |
48 | | - Hlr += 2.0, "Z", (2, 2), "Z", (3, 2) |
49 | | - Hlr += -1.0, "Z", (1, 2), "Z", (3, 1) |
50 | | - |
51 | | - # root_vertex = (1, 2) |
52 | | - # println(leaf_vertices(is)) |
53 | | - |
54 | | - @testset "Svd approach" for root_vertex in leaf_vertices(is) |
55 | | - # get TTN Hamiltonian directly |
56 | | - Hsvd = ttn(H, is; root_vertex, cutoff = 1.0e-10) |
57 | | - # get corresponding MPO Hamiltonian |
58 | | - Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], H), sites) |
59 | | - # compare resulting dense Hamiltonians |
60 | | - @disable_warn_order begin |
61 | | - Tttno = prod(Hline) |
62 | | - Tmpo = contract(Hsvd) |
63 | | - end |
64 | | - @test Tttno ≈ Tmpo rtol = 1.0e-6 |
65 | | - |
66 | | - Hsvd_lr = ttn(Hlr, is; root_vertex, cutoff = 1.0e-10) |
67 | | - Hline_lr = ITensorMPS.MPO(replace_vertices(v -> vmap[v], Hlr), sites) |
68 | | - @disable_warn_order begin |
69 | | - Tttno_lr = prod(Hline_lr) |
70 | | - Tmpo_lr = contract(Hsvd_lr) |
71 | | - end |
72 | | - @test Tttno_lr ≈ Tmpo_lr rtol = 1.0e-6 |
73 | | - end |
74 | | - if auto_fermion_enabled |
75 | | - ITensors.enable_auto_fermion() |
76 | | - end |
77 | | - end |
78 | | - |
79 | 9 | @testset "Multiple onsite terms (regression test for issue #62)" begin |
80 | 10 | auto_fermion_enabled = ITensors.using_auto_fermion() |
81 | 11 | if !auto_fermion_enabled |
|
98 | 28 | ITensors.enable_auto_fermion() |
99 | 29 | end |
100 | 30 | end |
101 | | - |
102 | | - @testset "OpSum to TTN QN" begin |
103 | | - # small comb tree |
104 | | - tooth_lengths = fill(2, 3) |
105 | | - c = named_comb_tree(tooth_lengths) |
106 | | - is = siteinds("S=1/2", c; conserve_qns = true) |
107 | | - is_noqns = copy(is) |
108 | | - for v in vertices(is) |
109 | | - is_noqns[v] = removeqns(is_noqns[v]) |
110 | | - end |
111 | | - |
112 | | - # linearized version |
113 | | - linear_order = [4, 1, 2, 5, 3, 6] |
114 | | - vmap = Dictionary(collect(vertices(is))[linear_order], eachindex(linear_order)) |
115 | | - sites = only.(collect(vertex_data(is)))[linear_order] |
116 | | - |
117 | | - # test with next-to-nearest-neighbor Ising Hamiltonian |
118 | | - J1 = -1 |
119 | | - J2 = 2 |
120 | | - h = 0.5 |
121 | | - H = ModelHamiltonians.heisenberg(c; J1 = J1, J2 = J2, h = h) |
122 | | - # add combination of longer range interactions |
123 | | - Hlr = copy(H) |
124 | | - Hlr += 5, "Z", (1, 2), "Z", (2, 2) #, "Z", (3,2) |
125 | | - Hlr += -4, "Z", (1, 1), "Z", (2, 2) |
126 | | - Hlr += 2.0, "Z", (2, 2), "Z", (3, 2) |
127 | | - Hlr += -1.0, "Z", (1, 2), "Z", (3, 1) |
128 | | - |
129 | | - # root_vertex = (1, 2) |
130 | | - # println(leaf_vertices(is)) |
131 | | - |
132 | | - @testset "Svd approach" for root_vertex in leaf_vertices(is) |
133 | | - # get TTN Hamiltonian directly |
134 | | - Hsvd = ttn(H, is; root_vertex, cutoff = 1.0e-10) |
135 | | - # get corresponding MPO Hamiltonian |
136 | | - Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], H), sites) |
137 | | - # compare resulting sparse Hamiltonians |
138 | | - |
139 | | - @disable_warn_order begin |
140 | | - Tmpo = prod(Hline) |
141 | | - Tttno = contract(Hsvd) |
142 | | - end |
143 | | - @test Tttno ≈ Tmpo rtol = 1.0e-6 |
144 | | - |
145 | | - Hsvd_lr = ttn(Hlr, is; root_vertex, cutoff = 1.0e-10) |
146 | | - Hline_lr = ITensorMPS.MPO(replace_vertices(v -> vmap[v], Hlr), sites) |
147 | | - @disable_warn_order begin |
148 | | - Tttno_lr = prod(Hline_lr) |
149 | | - Tmpo_lr = contract(Hsvd_lr) |
150 | | - end |
151 | | - @test Tttno_lr ≈ Tmpo_lr rtol = 1.0e-6 |
152 | | - end |
153 | | - end |
154 | | - |
155 | | - @testset "OpSum to TTN Fermions" begin |
156 | | - # small comb tree |
157 | | - auto_fermion_enabled = ITensors.using_auto_fermion() |
158 | | - if !auto_fermion_enabled |
159 | | - ITensors.enable_auto_fermion() |
160 | | - end |
161 | | - tooth_lengths = fill(2, 3) |
162 | | - c = named_comb_tree(tooth_lengths) |
163 | | - is = siteinds("Fermion", c; conserve_nf = true) |
164 | | - |
165 | | - # test with next-nearest neighbor tight-binding model |
166 | | - t = 1.0 |
167 | | - tp = 0.4 |
168 | | - U = 0.0 |
169 | | - h = 0.5 |
170 | | - H = ModelHamiltonians.tight_binding(c; t, tp, h) |
171 | | - |
172 | | - # add combination of longer range interactions |
173 | | - Hlr = copy(H) |
174 | | - |
175 | | - @testset "Svd approach" for root_vertex in leaf_vertices(is) |
176 | | - # get TTN Hamiltonian directly |
177 | | - Hsvd = ttn(H, is; root_vertex, cutoff = 1.0e-10) |
178 | | - # get corresponding MPO Hamiltonian |
179 | | - sites = [only(is[v]) for v in reverse(post_order_dfs_vertices(c, root_vertex))] |
180 | | - vmap = Dictionary( |
181 | | - reverse(post_order_dfs_vertices(c, root_vertex)), |
182 | | - 1:length(sites) |
183 | | - ) |
184 | | - Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], H), sites) |
185 | | - @disable_warn_order begin |
186 | | - Tmpo = prod(Hline) |
187 | | - Tttno = contract(Hsvd) |
188 | | - end |
189 | | - |
190 | | - # verify that the norm isn't 0 and thus the same (which would indicate a problem with the autofermion system |
191 | | - @test norm(Tmpo) > 0 |
192 | | - @test norm(Tttno) > 0 |
193 | | - @test norm(Tmpo) ≈ norm(Tttno) rtol = 1.0e-6 |
194 | | - |
195 | | - # TODO: fix comparison for fermionic tensors |
196 | | - @test_broken Tmpo ≈ Tttno |
197 | | - # In the meantime: matricize tensors and convert to dense Matrix to compare element by element |
198 | | - dTmm = to_matrix(Tmpo) |
199 | | - dTtm = to_matrix(Tttno) |
200 | | - @test any(>(1.0e-14), dTmm - dTtm) |
201 | | - end |
202 | | - if !auto_fermion_enabled |
203 | | - ITensors.disable_auto_fermion() |
204 | | - end |
205 | | - end |
206 | | - |
207 | | - @testset "OpSum to TTN QN missing" begin |
208 | | - # small comb tree |
209 | | - tooth_lengths = fill(2, 3) |
210 | | - c = named_comb_tree(tooth_lengths) |
211 | | - c2 = copy(c) |
212 | | - ## add an internal vertex into the comb graph c2 |
213 | | - add_vertex!(c2, (-1, 1)) |
214 | | - add_edge!(c2, (-1, 1) => (2, 2)) |
215 | | - add_edge!(c2, (-1, 1) => (3, 1)) |
216 | | - add_edge!(c2, (-1, 1) => (2, 1)) |
217 | | - rem_edge!(c2, (2, 1) => (2, 2)) |
218 | | - rem_edge!(c2, (2, 1) => (3, 1)) |
219 | | - |
220 | | - is = siteinds("S=1/2", c; conserve_qns = true) |
221 | | - is_missing_site = siteinds("S=1/2", c2; conserve_qns = true) |
222 | | - is_missing_site[(-1, 1)] = Vector{Index}[] |
223 | | - |
224 | | - # linearized version |
225 | | - linear_order = [4, 1, 2, 5, 3, 6] |
226 | | - vmap = Dictionary(collect(vertices(is))[linear_order], eachindex(linear_order)) |
227 | | - sites = |
228 | | - only.(filter(d -> !isempty(d), collect(vertex_data(is_missing_site))))[linear_order] |
229 | | - |
230 | | - J1 = -1 |
231 | | - J2 = 2 |
232 | | - h = 0.5 |
233 | | - # connectivity of the Hamiltonian is that of the original comb graph |
234 | | - H = ModelHamiltonians.heisenberg(c; J1, J2, h) |
235 | | - |
236 | | - # add combination of longer range interactions |
237 | | - Hlr = copy(H) |
238 | | - Hlr += 5, "Z", (1, 2), "Z", (2, 2) |
239 | | - Hlr += -4, "Z", (1, 1), "Z", (2, 2) |
240 | | - Hlr += 2.0, "Z", (2, 2), "Z", (3, 2) |
241 | | - Hlr += -1.0, "Z", (1, 2), "Z", (3, 1) |
242 | | - |
243 | | - @testset "Svd approach" for root_vertex in leaf_vertices(is) |
244 | | - # get TTN Hamiltonian directly |
245 | | - Hsvd = ttn(H, is_missing_site; root_vertex, cutoff = 1.0e-10) |
246 | | - # get corresponding MPO Hamiltonian |
247 | | - Hline = ITensorMPS.MPO(replace_vertices(v -> vmap[v], H), sites) |
248 | | - |
249 | | - # compare resulting sparse Hamiltonians |
250 | | - @disable_warn_order begin |
251 | | - Tmpo = prod(Hline) |
252 | | - Tttno = contract(Hsvd) |
253 | | - end |
254 | | - @test Tttno ≈ Tmpo rtol = 1.0e-6 |
255 | | - |
256 | | - Hsvd_lr = ttn(Hlr, is_missing_site; root_vertex, cutoff = 1.0e-10) |
257 | | - Hline_lr = ITensorMPS.MPO(replace_vertices(v -> vmap[v], Hlr), sites) |
258 | | - @disable_warn_order begin |
259 | | - Tttno_lr = prod(Hline_lr) |
260 | | - Tmpo_lr = contract(Hsvd_lr) |
261 | | - end |
262 | | - @test Tttno_lr ≈ Tmpo_lr rtol = 1.0e-6 |
263 | | - end |
264 | | - end |
265 | 31 | end |
266 | 32 | end |
0 commit comments