diff --git a/Project.toml b/Project.toml index 40f228c..23393dd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensorMPOConstruction" uuid = "f1f99f3b-4b00-4d2f-a9c2-ce76efdc8f31" -authors = ["Ben Corbett and contributors"] version = "0.1.0" +authors = ["Ben Corbett and contributors"] [deps] FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e" @@ -10,6 +10,7 @@ ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" @@ -17,6 +18,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" ITensorMPS = "0.1, 0.2, 0.3" ITensors = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9" LinearAlgebra = "1.9, 1.10" +Random = "1" SparseArrays = "1.9, 1.10" TimerOutputs = "0.5" julia = "1.10" diff --git a/src/MPOConstruction.jl b/src/MPOConstruction.jl index 8dae50b..192980f 100644 --- a/src/MPOConstruction.jl +++ b/src/MPOConstruction.jl @@ -89,7 +89,7 @@ function resume_svd_MPO( llinks[n + 1], prime(sites[n]), dag(sites[n]); - tol=1e-10, + tol=0.0, checkflux=false, ) else @@ -113,7 +113,7 @@ function resume_svd_MPO( llinks[n + 1], prime(sites[n]), dag(sites[n]); - tol=1e-10, + tol=0.0, checkflux=false, ) end diff --git a/test/data/h1_mponew_bad_case_terms.txt b/test/data/h1_mponew_bad_case_terms.txt new file mode 100644 index 0000000..4cf2aa4 --- /dev/null +++ b/test/data/h1_mponew_bad_case_terms.txt @@ -0,0 +1,253 @@ +# N=100 +1 25 37931.36296516507 +2 26 117103.84155871464 +3 27 117103.84155871464 +4 28 117103.84155871464 +5 25 -29545.40960407293 +6 26 -53236.781466396234 +7 27 -53236.781466396234 +8 28 -53236.781466396234 +9 25 23838.43229337311 +10 26 34835.93267151127 +11 27 34835.93267151127 +12 28 34835.93267151127 +13 25 -19263.286003783112 +14 26 -25226.711735954024 +15 27 -25226.711735954024 +16 28 -25226.711735954024 +17 25 22599.0145640706 +18 26 26061.170772204554 +19 27 26061.170772204554 +20 28 26061.170772204554 +21 25 -115532.52988556704 +21 29 19157.360826136017 +21 33 -3714.19421163894 +21 37 177.33414769923584 +21 41 -38.211839548830376 +21 45 139.23602772855048 +21 49 -88.15570397269616 +21 53 26.060553298381105 +21 57 1.331354445228164 +21 61 -4.498057184472989 +21 65 1.8189826521214016 +21 69 -0.3472052816509936 +21 73 0.030522533664478185 +21 77 -0.0014393368948888122 +21 81 -0.00012809052369549756 +21 85 0.00012429419120491152 +21 89 -1.9838069548747782e-5 +21 93 3.0918509292993344e-6 +21 97 3.866050749666841e-7 +22 26 -120349.02154222567 +22 30 20864.082202176112 +22 34 -4277.885131052926 +22 38 428.16342372165354 +22 42 -221.33807947417543 +22 46 251.0088012710466 +22 50 -141.43189174348248 +22 54 44.36481268487012 +22 58 -2.4989404631691334 +22 62 -4.29944651220679 +22 66 1.9628425431439491 +22 70 -0.39952395207949065 +22 74 0.04143748276920632 +22 78 -0.0036420864531109642 +22 82 0.0003376175177695186 +22 86 5.1606272815145875e-5 +22 90 -1.221093678762559e-5 +22 94 2.9991011840707995e-6 +22 98 4.72312404103241e-7 +23 27 -120349.02154222567 +23 31 20864.082202176112 +23 35 -4277.885131052926 +23 39 428.16342372165354 +23 43 -221.33807947417543 +23 47 251.0088012710466 +23 51 -141.43189174348248 +23 55 44.36481268487012 +23 59 -2.4989404631691334 +23 63 -4.29944651220679 +23 67 1.9628425431439491 +23 71 -0.39952395207949065 +23 75 0.04143748276920632 +23 79 -0.0036420864531109642 +23 83 0.0003376175177695186 +23 87 5.1606272815145875e-5 +23 91 -1.221093678762559e-5 +23 95 2.9991011840707995e-6 +23 99 4.72312404103241e-7 +24 28 -120349.02154222567 +24 32 20864.082202176112 +24 36 -4277.885131052926 +24 40 428.16342372165354 +24 44 -221.33807947417543 +24 48 251.0088012710466 +24 52 -141.43189174348248 +24 56 44.36481268487012 +24 60 -2.4989404631691334 +24 64 -4.29944651220679 +24 68 1.9628425431439491 +24 72 -0.39952395207949065 +24 76 0.04143748276920632 +24 80 -0.0036420864531109642 +24 84 0.0003376175177695186 +24 88 5.1606272815145875e-5 +24 92 -1.221093678762559e-5 +24 96 2.9991011840707995e-6 +24 100 4.72312404103241e-7 +25 1 37931.36296516507 +25 5 -29545.40960407293 +25 9 23838.43229337311 +25 13 -19263.286003783112 +25 17 22599.0145640706 +25 21 -115532.52988556704 +25 29 -56020.29866930369 +25 33 8878.808057849386 +25 37 -1388.426309729835 +25 41 306.34151560392195 +25 45 -261.1726800386604 +25 49 152.16482701959447 +25 53 -42.999742515087156 +25 57 -4.620607431348768 +26 2 117103.84155871464 +26 6 -53236.781466396234 +26 10 34835.93267151127 +26 14 -25226.711735954024 +26 18 26061.170772204554 +26 22 -120349.02154222567 +26 30 -57337.07910840612 +26 34 9248.937132581366 +26 38 -1497.2685821741113 +26 42 398.36880172003816 +26 46 -327.21618079577706 +26 50 185.51263652139676 +26 54 -54.589572263184735 +26 58 -2.299780204094004 +27 3 117103.84155871464 +27 7 -53236.781466396234 +27 11 34835.93267151127 +27 15 -25226.711735954024 +27 19 26061.170772204554 +27 23 -120349.02154222567 +27 31 -57337.07910840612 +27 35 9248.937132581366 +27 39 -1497.2685821741113 +27 43 398.36880172003816 +27 47 -327.21618079577706 +27 51 185.51263652139676 +27 55 -54.589572263184735 +28 4 117103.84155871464 +28 8 -53236.781466396234 +28 12 34835.93267151127 +28 16 -25226.711735954024 +28 20 26061.170772204554 +28 24 -120349.02154222567 +28 32 -57337.07910840612 +28 36 9248.937132581366 +28 40 -1497.2685821741113 +28 44 398.36880172003816 +28 48 -327.21618079577706 +28 52 185.51263652139676 +28 56 -54.589572263184735 +29 21 19157.360826136017 +29 25 -56020.29866930369 +30 22 20864.082202176112 +30 26 -57337.07910840612 +31 23 20864.082202176112 +31 27 -57337.07910840612 +32 24 20864.082202176112 +32 28 -57337.07910840612 +33 21 -3714.19421163894 +33 25 8878.808057849386 +34 22 -4277.885131052926 +34 26 9248.937132581366 +35 23 -4277.885131052926 +35 27 9248.937132581366 +36 24 -4277.885131052926 +36 28 9248.937132581366 +37 21 177.33414769923584 +37 25 -1388.426309729835 +38 22 428.16342372165354 +38 26 -1497.2685821741113 +39 23 428.16342372165354 +39 27 -1497.2685821741113 +40 24 428.16342372165354 +40 28 -1497.2685821741113 +41 21 -38.211839548830376 +41 25 306.34151560392195 +42 22 -221.33807947417543 +42 26 398.36880172003816 +43 23 -221.33807947417543 +43 27 398.36880172003816 +44 24 -221.33807947417543 +44 28 398.36880172003816 +45 21 139.23602772855048 +45 25 -261.1726800386604 +46 22 251.0088012710466 +46 26 -327.21618079577706 +47 23 251.0088012710466 +47 27 -327.21618079577706 +48 24 251.0088012710466 +48 28 -327.21618079577706 +49 21 -88.15570397269616 +49 25 152.16482701959447 +50 22 -141.43189174348248 +50 26 185.51263652139676 +51 23 -141.43189174348248 +51 27 185.51263652139676 +52 24 -141.43189174348248 +52 28 185.51263652139676 +53 21 26.060553298381105 +53 25 -42.999742515087156 +54 22 44.36481268487012 +54 26 -54.589572263184735 +55 23 44.36481268487012 +55 27 -54.589572263184735 +56 24 44.36481268487012 +56 28 -54.589572263184735 +57 21 1.331354445228164 +57 25 -4.620607431348768 +58 22 -2.4989404631691334 +58 26 -2.299780204094004 +59 23 -2.4989404631691334 +60 24 -2.4989404631691334 +61 21 -4.498057184472989 +62 22 -4.29944651220679 +63 23 -4.29944651220679 +64 24 -4.29944651220679 +65 21 1.8189826521214016 +66 22 1.9628425431439491 +67 23 1.9628425431439491 +68 24 1.9628425431439491 +69 21 -0.3472052816509936 +70 22 -0.39952395207949065 +71 23 -0.39952395207949065 +72 24 -0.39952395207949065 +73 21 0.030522533664478185 +74 22 0.04143748276920632 +75 23 0.04143748276920632 +76 24 0.04143748276920632 +77 21 -0.0014393368948888122 +78 22 -0.0036420864531109642 +79 23 -0.0036420864531109642 +80 24 -0.0036420864531109642 +81 21 -0.00012809052369549756 +82 22 0.0003376175177695186 +83 23 0.0003376175177695186 +84 24 0.0003376175177695186 +85 21 0.00012429419120491152 +86 22 5.1606272815145875e-5 +87 23 5.1606272815145875e-5 +88 24 5.1606272815145875e-5 +89 21 -1.9838069548747782e-5 +90 22 -1.221093678762559e-5 +91 23 -1.221093678762559e-5 +92 24 -1.221093678762559e-5 +93 21 3.0918509292993344e-6 +94 22 2.9991011840707995e-6 +95 23 2.9991011840707995e-6 +96 24 2.9991011840707995e-6 +97 21 3.866050749666841e-7 +98 22 4.72312404103241e-7 +99 23 4.72312404103241e-7 +100 24 4.72312404103241e-7 diff --git a/test/mponew_h1_accuracy.jl b/test/mponew_h1_accuracy.jl new file mode 100644 index 0000000..1e4ccb9 --- /dev/null +++ b/test/mponew_h1_accuracy.jl @@ -0,0 +1,141 @@ +using ITensorMPOConstruction +using ITensors +using ITensorMPS: siteinds, productMPS, MPO, OpSum, inner, apply, normalize! +using LinearAlgebra +using Random +using Test + +function read_h1_terms(path::AbstractString) + isfile(path) || error("terms file not found: $path") + N = 0 + terms = Tuple{Int,Int,Float64}[] + + for ln in eachline(path) + s = strip(ln) + isempty(s) && continue + + if startswith(s, "#") + if startswith(s, "# N=") + N = parse(Int, strip(s[5:end])) + end + continue + end + + f = split(s) + length(f) == 3 || error("invalid terms line: $ln") + p = parse(Int, f[1]) + q = parse(Int, f[2]) + h = parse(Float64, f[3]) + push!(terms, (p, q, h)) + N = max(N, p, q) + end + + isempty(terms) && error("no terms found in $path") + return N, terms +end + +function dense_h1(N::Int, terms) + H = zeros(Float64, N, N) + for (p, q, h) in terms + H[p, q] += h + end + return H +end + +function sample_sites(n::Int, k::Int) + k = max(1, min(k, n)) + idx = sort(unique(round.(Int, collect(range(1, n, length=k))))) + while length(idx) < k + for i in 1:n + (i in idx) && continue + push!(idx, i) + length(idx) == k && break + end + end + sort!(idx) + return idx +end + +function make_test_orbitals(N::Int; n_basis::Int=8, n_random::Int=4, seed::Int=1234) + tests = Vector{Vector{Float64}}() + + for p in sample_sites(N, n_basis) + phi = zeros(Float64, N) + phi[p] = 1.0 + push!(tests, phi) + end + + rng = MersenneTwister(seed) + for _ in 1:n_random + phi = randn(rng, N) + nrm = norm(phi) + nrm > 0 || (phi[1] = 1.0; nrm = 1.0) + phi ./= nrm + push!(tests, phi) + end + + return tests +end + +function onep_state_fermion(sites, phi::AbstractVector{<:Real}; tol::Real=1e-14) + N = length(sites) + length(phi) == N || error("orbital length mismatch") + + phi_n = Float64.(phi) + phi_n ./= norm(phi_n) + + p0 = argmax(abs.(phi_n)) + st = fill("Emp", N) + st[p0] = "Occ" + psi0 = productMPS(sites, st) + + os = OpSum() + for i in 1:N + c = phi_n[i] + abs(c) <= tol && continue + if i == p0 + os += c, "N", p0 + else + os += c, "Cdag", i, "C", p0 + end + end + + psi = apply(MPO(os, sites), psi0; alg="naive", cutoff=0.0, maxdim=typemax(Int)) + normalize!(psi) + return psi +end + +function build_h1_opsum(terms) + os = OpSum() + for (p, q, h) in terms + if p == q + os += h, "N", p + else + os += h, "Cdag", p, "C", q + end + end + return os +end + +@testset "MPO_new H1 accuracy regression" begin + data_path = joinpath(@__DIR__, "data", "h1_mponew_bad_case_terms.txt") + N, terms = read_h1_terms(data_path) + + H = dense_h1(N, terms) + @test maximum(abs, H .- transpose(H)) <= 1e-12 + + sites = siteinds("Fermion", N; conserve_qns=true, conserve_nf=true) + H_new = MPO_new(build_h1_opsum(terms), sites; tol=1.0) + + tests = make_test_orbitals(N; n_basis=8, n_random=4, seed=1234) + max_delta = 0.0 + for phi in tests + psi = onep_state_fermion(sites, phi) + E_dense = real(dot(phi, H * phi)) + E_new = real(inner(psi', H_new, psi)) + max_delta = max(max_delta, abs(E_new - E_dense)) + end + + # Before fixing the hard-coded assembly cutoff this was ~3e-7 for this case. + @test max_delta <= 1e-9 +end diff --git a/test/runtests.jl b/test/runtests.jl index f7042fc..b8ddced 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,5 @@ include("large-graph-tests.jl") include("ops-tests.jl") include("test-MPOConstruction.jl") +include("mponew_h1_accuracy.jl") include("../examples/haldane-shastry.jl")