diff --git a/Project.toml b/Project.toml index b649d2cb..a989b9b4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ITensorNetworks" uuid = "2919e153-833c-4bdc-8836-1ea460a35fc7" authors = ["Matthew Fishman , Joseph Tindall and contributors"] -version = "0.13.15" +version = "0.13.16" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" diff --git a/src/lib/ITensorsExtensions/src/itensor.jl b/src/lib/ITensorsExtensions/src/itensor.jl index 8c834830..cc03865f 100644 --- a/src/lib/ITensorsExtensions/src/itensor.jl +++ b/src/lib/ITensorsExtensions/src/itensor.jl @@ -1,16 +1,20 @@ using LinearAlgebra: LinearAlgebra, eigen, pinv using ITensors: + ITensors, ITensor, Index, commonind, dag, + dir, hasqns, + indices, inds, isdiag, itensor, map_diag, noncommonind, noprime, + permute, replaceind, replaceinds, sim, @@ -58,18 +62,57 @@ function eigendecomp(A::ITensor, linds, rinds; ishermitian=false, kwargs...) D, U = eigen(A, linds, rinds; ishermitian, kwargs...) ul, ur = noncommonind(D, U), commonind(D, U) Ul = replaceinds(U, vcat(rinds, ur), vcat(linds, ul)) - return Ul, D, dag(U) end -function map_eigvals(f::Function, A::ITensor, inds...; ishermitian=false, kwargs...) +function make_bosonic(A::ITensor, Linds, Rinds) + # Bring indices into i',j',..,dag(j),dag(i) + # ordering with Out indices coming before In indices + # Resulting tensor acts like a normal matrix (no extra signs + # when taking powers A^n) + Linds = indices(Linds) + Rinds = indices(Rinds) + if all(j->dir(j)==ITensors.Out, Linds) && all(j->dir(j)==ITensors.In, Rinds) + ordered_inds = [Linds..., reverse(Rinds)...] + elseif all(j->dir(j)==ITensors.Out, Rinds) && all(j->dir(j)==ITensors.In, Linds) + ordered_inds = [Rinds..., reverse(Linds)...] + else + error( + "For fermionic exp, Linds and Rinds must have same directions within each set. Got dir.(Linds)=", + dir.(Linds), + ", dir.(Rinds)=", + dir.(Rinds), + ) + end + # permuted A^n will be sign free, ok to temporarily disable fermion system + return permute(A, ordered_inds) +end + +function map_eigvals(f::Function, A::ITensor, Linds, Rinds; kws...) + # + fermionic_itensor = + ITensors.using_auto_fermion() && ITensors.has_fermionic_subspaces(inds(A)) + if fermionic_itensor + A = make_bosonic(A::ITensor, Linds, Rinds) + ordered_inds = inds(A) + ITensors.disable_auto_fermion() + end + if isdiag(A) - return map_diag(f, A) + mapped_A = map_diag(f, A) + else + Ul, D, Ur = eigendecomp(A, Linds, Rinds; kws...) + mapped_A = Ul * map_diag(f, D) * Ur end - Ul, D, Ur = eigendecomp(A, inds...; ishermitian, kwargs...) + # + if fermionic_itensor + # Ensure indices in "matrix" form before re-enabling fermion system + mapped_A = permute(mapped_A, ordered_inds) + ITensors.enable_auto_fermion() + end - return Ul * map_diag(f, D) * Ur + return mapped_A end # Analagous to `denseblocks`. diff --git a/test/test_itensorsextensions.jl b/test/test_itensorsextensions.jl index cb4b5ed0..39ab8005 100644 --- a/test/test_itensorsextensions.jl +++ b/test/test_itensorsextensions.jl @@ -4,17 +4,22 @@ using ITensors: ITensor, Index, QN, + apply, dag, delta, inds, + mapprime, noprime, + norm, op, + permute, prime, random_itensor, replaceind, replaceinds, - sim -using ITensorNetworks.ITensorsExtensions: map_eigvals + sim, + swapprime +using ITensorNetworks.ITensorsExtensions: eigendecomp, map_eigvals using StableRNGs: StableRNG using Test: @test, @testset @testset "ITensorsExtensions" begin @@ -55,9 +60,9 @@ using Test: @test, @testset rng = StableRNG(1234) A = random_itensor(rng, elt, i, j) P = A * prime(dag(A), i) - sqrtP = map_eigvals(sqrt, P, i, i'; ishermitian=true) - inv_P = dag(map_eigvals(inv, P, i, i'; ishermitian=true)) - inv_sqrtP = dag(map_eigvals(inv ∘ sqrt, P, i, i'; ishermitian=true)) + sqrtP = map_eigvals(sqrt, P, [i], [i']; ishermitian=true) + inv_P = dag(map_eigvals(inv, P, [i], [i']; ishermitian=true)) + inv_sqrtP = dag(map_eigvals(inv ∘ sqrt, P, [i], [i']; ishermitian=true)) new_ind = noprime(sim(i')) sqrtPdag = replaceind(dag(sqrtP), i', new_ind) @@ -72,5 +77,66 @@ using Test: @test, @testset I = replaceind(inv_sqrtP * sqrtP, new_ind, i) @test I ≈ op("I", i) end + + @testset "Fermionic eigendecomp" begin + ITensors.enable_auto_fermion() + s1 = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "Site,Fermion,n=1") + s2 = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "Site,Fermion,n=2") + + # Make a random Hermitian matrix-like 4th order ITensor + T = random_itensor(s1', s2', dag(s2), dag(s1)) + T = apply(T, swapprime(dag(T), 0=>1)) + @test T ≈ swapprime(dag(T), 0=>1) # check Hermitian + + Ul, D, Ur = eigendecomp(T, [s1', s2'], [dag(s1), dag(s2)]; ishermitian=true) + + @test Ul*D*Ur ≈ T + ITensors.disable_auto_fermion() + end + + @testset "Fermionic map eigvals tests" begin + ITensors.enable_auto_fermion() + s1 = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "Site,Fermion,n=1") + s2 = Index([QN("Nf", 0, -1)=>2, QN("Nf", 1, -1)=>2], "Site,Fermion,n=2") + + # Make a random Hermitian matrix ITensor + M = random_itensor(s1', dag(s1)) + #M = mapprime(prime(M)*swapprime(dag(M),0=>1),2=>1) + M = apply(M, swapprime(dag(M), 0=>1)) + + # Make a random Hermitian matrix-like 4th order ITensor + T = random_itensor(s1', s2', dag(s2), dag(s1)) + T = apply(T, swapprime(dag(T), 0=>1)) + + # Matrix test + sqrtM = map_eigvals(sqrt, M, [s1'], [dag(s1)]; ishermitian=true) + @test M ≈ apply(sqrtM, sqrtM) + + ## Tensor test + sqrtT = map_eigvals(sqrt, T, [s1', s2'], [dag(s1), dag(s2)]; ishermitian=true) + @test T ≈ apply(sqrtT, sqrtT) + + # Permute and test again + T = permute(T, dag(s2), s2', dag(s1), s1') + sqrtT = map_eigvals(sqrt, T, [s1', s2'], [dag(s1), dag(s2)]; ishermitian=true) + @test T ≈ apply(sqrtT, sqrtT) + + ## Explicitly passing indices in different, valid orders + sqrtT = map_eigvals(sqrt, T, [s2', s1'], [dag(s2), dag(s1)]; ishermitian=true) + @test T ≈ apply(sqrtT, sqrtT) + sqrtT = map_eigvals(sqrt, T, [dag(s2), dag(s1)], [s2', s1'], ; ishermitian=true) + @test T ≈ apply(sqrtT, sqrtT) + sqrtT = map_eigvals(sqrt, T, [dag(s1), dag(s2)], [s1', s2'], ; ishermitian=true) + @test T ≈ apply(sqrtT, sqrtT) + + # Test bosonic index case while fermion system is enabled + b = Index([QN("Nb", 0)=>2, QN("Nb", 1)=>2]) + T = random_itensor(b', dag(b)) + T = apply(T, swapprime(dag(T), 0=>1)) + sqrtT = map_eigvals(sqrt, T, [b'], [dag(b)]; ishermitian=true) + @test T ≈ apply(sqrtT, sqrtT) + + ITensors.disable_auto_fermion() + end end end