From fc2519479794db92d7e3cf4dcebcc20ef1228298 Mon Sep 17 00:00:00 2001 From: Miles Date: Tue, 29 Mar 2022 23:12:12 -0400 Subject: [PATCH 1/3] Add time-dependent Green functions --- ITensorGaussianMPS/src/ITensorGaussianMPS.jl | 3 + ITensorGaussianMPS/src/dynamics.jl | 88 ++++++++++++++++++++ ITensorGaussianMPS/test/dynamics.jl | 55 ++++++++++++ ITensorGaussianMPS/test/runtests.jl | 1 + 4 files changed, 147 insertions(+) create mode 100644 ITensorGaussianMPS/src/dynamics.jl create mode 100644 ITensorGaussianMPS/test/dynamics.jl diff --git a/ITensorGaussianMPS/src/ITensorGaussianMPS.jl b/ITensorGaussianMPS/src/ITensorGaussianMPS.jl index 7669a5645e..26b9b9f471 100644 --- a/ITensorGaussianMPS/src/ITensorGaussianMPS.jl +++ b/ITensorGaussianMPS/src/ITensorGaussianMPS.jl @@ -10,6 +10,9 @@ import LinearAlgebra: Givens export slater_determinant_to_mps, slater_determinant_to_gmps, hopping_hamiltonian, slater_determinant_matrix +export G_R,G_G,G_L + include("gmps.jl") +include("dynamics.jl") end diff --git a/ITensorGaussianMPS/src/dynamics.jl b/ITensorGaussianMPS/src/dynamics.jl new file mode 100644 index 0000000000..a59d1c0179 --- /dev/null +++ b/ITensorGaussianMPS/src/dynamics.jl @@ -0,0 +1,88 @@ + + +function _compute_correlator(N::Int, + f::Function; + kwargs...) + sites = get(kwargs,:sites,1:N) + site_range = (sites isa AbstractRange) ? sites : collect(sites) + Ns = length(site_range) + + C = zeros(ComplexF64,Ns,Ns) + for (ni,i) in enumerate(site_range), (nj,j) in enumerate(site_range) + C[ni,nj] += f(i,j) + end + + if sites isa Number + return C[1,1] + end + return C +end + +""" +Compute the retarded single-particle Green function +G_R(t)_ij = -i θ(t) <ϕ|{cᵢ(t), c†ⱼ(0)}|ϕ> +where ϕ is a Slater determinant +""" +function G_R(t::Number, + ϕ, + ϵ::Vector{Float64}; + kwargs...) + N = length(ϵ) + @assert size(ϕ) == (N,N) + function compute_GR(i,j) + gr = 0.0im + for n = 1:N + gr += -im*ϕ[i,n]*conj(ϕ[j,n])*exp(-im*ϵ[n]*t) + end + return gr + end + return _compute_correlator(N,compute_GR; kwargs...) +end + +""" +Compute the greater single-particle Green function +G_G(t)_ij = G>(t)_ij = -i <ϕ|cᵢ(t), c†ⱼ(0)|ϕ> +where ϕ is a Slater determinant +""" +function G_G(t::Number, + ϕ, + ϵ::Vector{Float64}; + kwargs...) + N = length(ϵ) + @assert size(ϕ) == (N,N) + Nneg = count(en->(en<0),ϵ) + Npart = get(kwargs,:Npart,Nneg) + + function compute_GG(i,j) + gg = 0.0im + for n = Npart+1:N + gg += -im*ϕ[i,n]*conj(ϕ[j,n])*exp(-im*ϵ[n]*t) + end + return gg + end + return _compute_correlator(N,compute_GG; kwargs...) +end + +""" +Compute the lesser single-particle Green function +G_L(t)_ij = G<(t)_ij = +i <ϕ|c†ᵢ(0), cⱼ(t)|ϕ> +where ϕ is a Slater determinant +""" +function G_L(t::Number, + ϕ, + ϵ::Vector{Float64}; + kwargs...) + N = length(ϵ) + @assert size(ϕ) == (N,N) + Nneg = count(en->(en<0),ϵ) + Npart = get(kwargs,:Npart,Nneg) + + function compute_GL(i,j) + gl = 0.0im + for n = 1:Npart + gl += +im*ϕ[i,n]*conj(ϕ[j,n])*exp(-im*ϵ[n]*t) + end + return gl + end + return _compute_correlator(N,compute_GL; kwargs...) +end diff --git a/ITensorGaussianMPS/test/dynamics.jl b/ITensorGaussianMPS/test/dynamics.jl new file mode 100644 index 0000000000..2b162c6dec --- /dev/null +++ b/ITensorGaussianMPS/test/dynamics.jl @@ -0,0 +1,55 @@ +using ITensorGaussianMPS +using ITensors +using LinearAlgebra +using Test + +function diagonalize_h(N) + hterms = OpSum() + for j=1:N-1 + hterms += -1,"Cdag",j,"C",j+1 + hterms += -1,"Cdag",j+1,"C",j + end + + h = hopping_hamiltonian(hterms) + ϵ, phi = eigen(h) + return ϵ,phi +end + +@testset "Green Function Properties" begin + N = 20 + ϵ,phi = diagonalize_h(N) + + t = 1.0 + + Gr = G_R(t,phi,ϵ) + Gl = G_L(t,phi,ϵ) + Gg = G_G(t,phi,ϵ) + + @test Gr ≈ (Gg-Gl) + + # Gr is -i at time zero: + @test G_R(0.0,phi,ϵ)[1,1] ≈ -1.0im +end + +@testset "Sites Keyword" begin + N = 20 + ϵ,phi = diagonalize_h(N) + + t = 1.0 + + G = G_R(t,phi,ϵ) + + G5_10 = G_R(t,phi,ϵ; sites=5:10) + @test size(G5_10) == (6,6) + @test G5_10 ≈ G[5:10,5:10] + + g = G_R(t,phi,ϵ; sites=7) + @test g isa Number + + # Non-contiguous case + Gnc = G_R(t,phi,ϵ; sites=[5,10,15]) + @test size(Gnc) == (3,3) +end + + +nothing diff --git a/ITensorGaussianMPS/test/runtests.jl b/ITensorGaussianMPS/test/runtests.jl index f3b54e441d..c9defcb015 100644 --- a/ITensorGaussianMPS/test/runtests.jl +++ b/ITensorGaussianMPS/test/runtests.jl @@ -5,4 +5,5 @@ using Test @testset "ITensorGaussianMPS.jl" begin include("gmps.jl") include("electron.jl") + include("dynamics.jl") end From d50a7d2741f592aa0e14d191c0ccf8ed1326a2f3 Mon Sep 17 00:00:00 2001 From: Miles Date: Tue, 29 Mar 2022 23:20:38 -0400 Subject: [PATCH 2/3] Formatting --- ITensorGaussianMPS/src/ITensorGaussianMPS.jl | 2 +- ITensorGaussianMPS/src/dynamics.jl | 68 ++++++++------------ ITensorGaussianMPS/test/dynamics.jl | 37 ++++++----- 3 files changed, 47 insertions(+), 60 deletions(-) diff --git a/ITensorGaussianMPS/src/ITensorGaussianMPS.jl b/ITensorGaussianMPS/src/ITensorGaussianMPS.jl index 26b9b9f471..07519fe9e0 100644 --- a/ITensorGaussianMPS/src/ITensorGaussianMPS.jl +++ b/ITensorGaussianMPS/src/ITensorGaussianMPS.jl @@ -10,7 +10,7 @@ import LinearAlgebra: Givens export slater_determinant_to_mps, slater_determinant_to_gmps, hopping_hamiltonian, slater_determinant_matrix -export G_R,G_G,G_L +export G_R, G_G, G_L include("gmps.jl") include("dynamics.jl") diff --git a/ITensorGaussianMPS/src/dynamics.jl b/ITensorGaussianMPS/src/dynamics.jl index a59d1c0179..73c884f55a 100644 --- a/ITensorGaussianMPS/src/dynamics.jl +++ b/ITensorGaussianMPS/src/dynamics.jl @@ -1,19 +1,16 @@ - -function _compute_correlator(N::Int, - f::Function; - kwargs...) - sites = get(kwargs,:sites,1:N) +function _compute_correlator(N::Int, f::Function; kwargs...) + sites = get(kwargs, :sites, 1:N) site_range = (sites isa AbstractRange) ? sites : collect(sites) Ns = length(site_range) - C = zeros(ComplexF64,Ns,Ns) - for (ni,i) in enumerate(site_range), (nj,j) in enumerate(site_range) - C[ni,nj] += f(i,j) + C = zeros(ComplexF64, Ns, Ns) + for (ni, i) in enumerate(site_range), (nj, j) in enumerate(site_range) + C[ni, nj] += f(i, j) end if sites isa Number - return C[1,1] + return C[1, 1] end return C end @@ -23,20 +20,17 @@ Compute the retarded single-particle Green function G_R(t)_ij = -i θ(t) <ϕ|{cᵢ(t), c†ⱼ(0)}|ϕ> where ϕ is a Slater determinant """ -function G_R(t::Number, - ϕ, - ϵ::Vector{Float64}; - kwargs...) +function G_R(t::Number, ϕ, ϵ::Vector{Float64}; kwargs...) N = length(ϵ) - @assert size(ϕ) == (N,N) - function compute_GR(i,j) + @assert size(ϕ) == (N, N) + function compute_GR(i, j) gr = 0.0im - for n = 1:N - gr += -im*ϕ[i,n]*conj(ϕ[j,n])*exp(-im*ϵ[n]*t) + for n in 1:N + gr += -im * ϕ[i, n] * conj(ϕ[j, n]) * exp(-im * ϵ[n] * t) end return gr end - return _compute_correlator(N,compute_GR; kwargs...) + return _compute_correlator(N, compute_GR; kwargs...) end """ @@ -44,23 +38,20 @@ Compute the greater single-particle Green function G_G(t)_ij = G>(t)_ij = -i <ϕ|cᵢ(t), c†ⱼ(0)|ϕ> where ϕ is a Slater determinant """ -function G_G(t::Number, - ϕ, - ϵ::Vector{Float64}; - kwargs...) +function G_G(t::Number, ϕ, ϵ::Vector{Float64}; kwargs...) N = length(ϵ) - @assert size(ϕ) == (N,N) - Nneg = count(en->(en<0),ϵ) - Npart = get(kwargs,:Npart,Nneg) + @assert size(ϕ) == (N, N) + Nneg = count(en -> (en < 0), ϵ) + Npart = get(kwargs, :Npart, Nneg) - function compute_GG(i,j) + function compute_GG(i, j) gg = 0.0im - for n = Npart+1:N - gg += -im*ϕ[i,n]*conj(ϕ[j,n])*exp(-im*ϵ[n]*t) + for n in (Npart + 1):N + gg += -im * ϕ[i, n] * conj(ϕ[j, n]) * exp(-im * ϵ[n] * t) end return gg end - return _compute_correlator(N,compute_GG; kwargs...) + return _compute_correlator(N, compute_GG; kwargs...) end """ @@ -68,21 +59,18 @@ Compute the lesser single-particle Green function G_L(t)_ij = G<(t)_ij = +i <ϕ|c†ᵢ(0), cⱼ(t)|ϕ> where ϕ is a Slater determinant """ -function G_L(t::Number, - ϕ, - ϵ::Vector{Float64}; - kwargs...) +function G_L(t::Number, ϕ, ϵ::Vector{Float64}; kwargs...) N = length(ϵ) - @assert size(ϕ) == (N,N) - Nneg = count(en->(en<0),ϵ) - Npart = get(kwargs,:Npart,Nneg) + @assert size(ϕ) == (N, N) + Nneg = count(en -> (en < 0), ϵ) + Npart = get(kwargs, :Npart, Nneg) - function compute_GL(i,j) + function compute_GL(i, j) gl = 0.0im - for n = 1:Npart - gl += +im*ϕ[i,n]*conj(ϕ[j,n])*exp(-im*ϵ[n]*t) + for n in 1:Npart + gl += +im * ϕ[i, n] * conj(ϕ[j, n]) * exp(-im * ϵ[n] * t) end return gl end - return _compute_correlator(N,compute_GL; kwargs...) + return _compute_correlator(N, compute_GL; kwargs...) end diff --git a/ITensorGaussianMPS/test/dynamics.jl b/ITensorGaussianMPS/test/dynamics.jl index 2b162c6dec..b18a4aa0b6 100644 --- a/ITensorGaussianMPS/test/dynamics.jl +++ b/ITensorGaussianMPS/test/dynamics.jl @@ -5,51 +5,50 @@ using Test function diagonalize_h(N) hterms = OpSum() - for j=1:N-1 - hterms += -1,"Cdag",j,"C",j+1 - hterms += -1,"Cdag",j+1,"C",j + for j in 1:(N - 1) + hterms += -1, "Cdag", j, "C", j + 1 + hterms += -1, "Cdag", j + 1, "C", j end h = hopping_hamiltonian(hterms) ϵ, phi = eigen(h) - return ϵ,phi + return ϵ, phi end @testset "Green Function Properties" begin N = 20 - ϵ,phi = diagonalize_h(N) + ϵ, phi = diagonalize_h(N) t = 1.0 - Gr = G_R(t,phi,ϵ) - Gl = G_L(t,phi,ϵ) - Gg = G_G(t,phi,ϵ) + Gr = G_R(t, phi, ϵ) + Gl = G_L(t, phi, ϵ) + Gg = G_G(t, phi, ϵ) - @test Gr ≈ (Gg-Gl) + @test Gr ≈ (Gg - Gl) # Gr is -i at time zero: - @test G_R(0.0,phi,ϵ)[1,1] ≈ -1.0im + @test G_R(0.0, phi, ϵ)[1, 1] ≈ -1.0im end @testset "Sites Keyword" begin N = 20 - ϵ,phi = diagonalize_h(N) + ϵ, phi = diagonalize_h(N) t = 1.0 - G = G_R(t,phi,ϵ) + G = G_R(t, phi, ϵ) - G5_10 = G_R(t,phi,ϵ; sites=5:10) - @test size(G5_10) == (6,6) - @test G5_10 ≈ G[5:10,5:10] + G5_10 = G_R(t, phi, ϵ; sites=5:10) + @test size(G5_10) == (6, 6) + @test G5_10 ≈ G[5:10, 5:10] - g = G_R(t,phi,ϵ; sites=7) + g = G_R(t, phi, ϵ; sites=7) @test g isa Number # Non-contiguous case - Gnc = G_R(t,phi,ϵ; sites=[5,10,15]) - @test size(Gnc) == (3,3) + Gnc = G_R(t, phi, ϵ; sites=[5, 10, 15]) + @test size(Gnc) == (3, 3) end - nothing From ba66adc5a20a07a1625e72106f536a25a38f2f58 Mon Sep 17 00:00:00 2001 From: Miles Date: Wed, 30 Mar 2022 22:42:38 -0400 Subject: [PATCH 3/3] Rename Green function functions and improve doc strings --- ITensorGaussianMPS/src/ITensorGaussianMPS.jl | 9 +++-- ITensorGaussianMPS/src/dynamics.jl | 35 +++++++++++++----- ITensorGaussianMPS/test/dynamics.jl | 38 ++++++++++++-------- 3 files changed, 56 insertions(+), 26 deletions(-) diff --git a/ITensorGaussianMPS/src/ITensorGaussianMPS.jl b/ITensorGaussianMPS/src/ITensorGaussianMPS.jl index 07519fe9e0..1d0745cdd4 100644 --- a/ITensorGaussianMPS/src/ITensorGaussianMPS.jl +++ b/ITensorGaussianMPS/src/ITensorGaussianMPS.jl @@ -8,9 +8,12 @@ using LinearAlgebra import LinearAlgebra: Givens export slater_determinant_to_mps, - slater_determinant_to_gmps, hopping_hamiltonian, slater_determinant_matrix - -export G_R, G_G, G_L + slater_determinant_to_gmps, + hopping_hamiltonian, + slater_determinant_matrix, + retarded_green_function, + lesser_green_function, + greater_green_function include("gmps.jl") include("dynamics.jl") diff --git a/ITensorGaussianMPS/src/dynamics.jl b/ITensorGaussianMPS/src/dynamics.jl index 73c884f55a..ee1793ace8 100644 --- a/ITensorGaussianMPS/src/dynamics.jl +++ b/ITensorGaussianMPS/src/dynamics.jl @@ -16,11 +16,14 @@ function _compute_correlator(N::Int, f::Function; kwargs...) end """ + retarded_green_function(t, ϕ, ϵ; kwargs...) + Compute the retarded single-particle Green function -G_R(t)_ij = -i θ(t) <ϕ|{cᵢ(t), c†ⱼ(0)}|ϕ> -where ϕ is a Slater determinant +GR(t)_ij = -i θ(t) <ϕ|{cᵢ(t), c†ⱼ(0)}|ϕ> +at time t where |ϕ> is a Slater determinant determined by +the matrix of orbitals ϕ, with energies ϵ. """ -function G_R(t::Number, ϕ, ϵ::Vector{Float64}; kwargs...) +function retarded_green_function(t::Number, ϕ, ϵ::Vector{Float64}; kwargs...) N = length(ϵ) @assert size(ϕ) == (N, N) function compute_GR(i, j) @@ -34,11 +37,18 @@ function G_R(t::Number, ϕ, ϵ::Vector{Float64}; kwargs...) end """ + greater_green_function(t, ϕ, ϵ; kwargs...) + Compute the greater single-particle Green function -G_G(t)_ij = G>(t)_ij = -i <ϕ|cᵢ(t), c†ⱼ(0)|ϕ> -where ϕ is a Slater determinant +G>(t)_ij = -i <ϕ|cᵢ(t) c†ⱼ(0)|ϕ> +at time t where |ϕ> is a Slater determinant determined by +the matrix of orbitals ϕ, with energies ϵ. + +The keyword argument `Npart` can be used to specify the number +of particles (number of occupied orbitals). Otherwise the number +of particles is determined by the number of negative energies in ϵ. """ -function G_G(t::Number, ϕ, ϵ::Vector{Float64}; kwargs...) +function greater_green_function(t::Number, ϕ, ϵ::Vector{Float64}; kwargs...) N = length(ϵ) @assert size(ϕ) == (N, N) Nneg = count(en -> (en < 0), ϵ) @@ -55,11 +65,18 @@ function G_G(t::Number, ϕ, ϵ::Vector{Float64}; kwargs...) end """ + lesser_green_function(t, ϕ, ϵ; kwargs...) + Compute the lesser single-particle Green function -G_L(t)_ij = G<(t)_ij = +i <ϕ|c†ᵢ(0), cⱼ(t)|ϕ> -where ϕ is a Slater determinant +G<(t)_ij = +i <ϕ|c†ᵢ(0) cⱼ(t)|ϕ> +at time t where |ϕ> is a Slater determinant determined by +the matrix of orbitals ϕ, with energies ϵ. + +The keyword argument `Npart` can be used to specify the number +of particles (number of occupied orbitals). Otherwise the number +of particles is determined by the number of negative energies in ϵ. """ -function G_L(t::Number, ϕ, ϵ::Vector{Float64}; kwargs...) +function lesser_green_function(t::Number, ϕ, ϵ::Vector{Float64}; kwargs...) N = length(ϵ) @assert size(ϕ) == (N, N) Nneg = count(en -> (en < 0), ϵ) diff --git a/ITensorGaussianMPS/test/dynamics.jl b/ITensorGaussianMPS/test/dynamics.jl index b18a4aa0b6..6364246c81 100644 --- a/ITensorGaussianMPS/test/dynamics.jl +++ b/ITensorGaussianMPS/test/dynamics.jl @@ -11,44 +11,54 @@ function diagonalize_h(N) end h = hopping_hamiltonian(hterms) - ϵ, phi = eigen(h) - return ϵ, phi + ϵ, ϕ = eigen(h) + return ϵ, ϕ end @testset "Green Function Properties" begin N = 20 - ϵ, phi = diagonalize_h(N) + ϵ, ϕ = diagonalize_h(N) t = 1.0 - Gr = G_R(t, phi, ϵ) - Gl = G_L(t, phi, ϵ) - Gg = G_G(t, phi, ϵ) + Gr = retarded_green_function(t, ϕ, ϵ) + Gl = lesser_green_function(t, ϕ, ϵ) + Gg = greater_green_function(t, ϕ, ϵ) @test Gr ≈ (Gg - Gl) - # Gr is -i at time zero: - @test G_R(0.0, phi, ϵ)[1, 1] ≈ -1.0im + # Gr is -im at time zero: + @test retarded_green_function(0.0, ϕ, ϵ)[1, 1] ≈ -im + + # G< at t=0 is im times the correlation matrix + Npart = div(N, 2) + Gl0 = lesser_green_function(0, ϕ, ϵ; Npart) + ϕ_occ = ϕ[:, 1:Npart] + C = ϕ_occ * ϕ_occ' + @test Gl0 ≈ im * C end @testset "Sites Keyword" begin N = 20 - ϵ, phi = diagonalize_h(N) + ϵ, ϕ = diagonalize_h(N) t = 1.0 - G = G_R(t, phi, ϵ) + G = retarded_green_function(t, ϕ, ϵ) - G5_10 = G_R(t, phi, ϵ; sites=5:10) + G5_10 = retarded_green_function(t, ϕ, ϵ; sites=5:10) @test size(G5_10) == (6, 6) @test G5_10 ≈ G[5:10, 5:10] - g = G_R(t, phi, ϵ; sites=7) - @test g isa Number + g7 = retarded_green_function(t, ϕ, ϵ; sites=7) + @test g7 isa Number # Non-contiguous case - Gnc = G_R(t, phi, ϵ; sites=[5, 10, 15]) + Gnc = retarded_green_function(t, ϕ, ϵ; sites=[5, 10, 15]) @test size(Gnc) == (3, 3) + @test Gnc[1, 1] ≈ G[5, 5] + @test Gnc[1, 2] ≈ G[5, 10] + @test Gnc[2, 3] ≈ G[10, 15] end nothing