diff --git a/ITensorGaussianMPS/src/ITensorGaussianMPS.jl b/ITensorGaussianMPS/src/ITensorGaussianMPS.jl
index 5885d83f60..ebf650247e 100644
--- a/ITensorGaussianMPS/src/ITensorGaussianMPS.jl
+++ b/ITensorGaussianMPS/src/ITensorGaussianMPS.jl
@@ -9,11 +9,15 @@ import LinearAlgebra: Givens
 
 export slater_determinant_to_mps,
   slater_determinant_to_gmps,
+  slater_determinant_to_gmera,
   hopping_hamiltonian,
   slater_determinant_matrix,
-  slater_determinant_to_gmera
+  retarded_green_function,
+  lesser_green_function,
+  greater_green_function
 
 include("gmps.jl")
 include("gmera.jl")
+include("dynamics.jl")
 
 end
diff --git a/ITensorGaussianMPS/src/dynamics.jl b/ITensorGaussianMPS/src/dynamics.jl
new file mode 100644
index 0000000000..ee1793ace8
--- /dev/null
+++ b/ITensorGaussianMPS/src/dynamics.jl
@@ -0,0 +1,93 @@
+
+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
+
+"""
+    retarded_green_function(t, ϕ, ϵ; kwargs...)
+
+Compute the retarded single-particle Green function
+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 retarded_green_function(t::Number, ϕ, ϵ::Vector{Float64}; kwargs...)
+  N = length(ϵ)
+  @assert size(ϕ) == (N, N)
+  function compute_GR(i, j)
+    gr = 0.0im
+    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...)
+end
+
+"""
+    greater_green_function(t, ϕ, ϵ; kwargs...)
+
+Compute the greater single-particle Green function
+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 greater_green_function(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 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...)
+end
+
+"""
+    lesser_green_function(t, ϕ, ϵ; kwargs...)
+
+Compute the lesser single-particle Green function
+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 lesser_green_function(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 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...)
+end
diff --git a/ITensorGaussianMPS/test/dynamics.jl b/ITensorGaussianMPS/test/dynamics.jl
new file mode 100644
index 0000000000..6364246c81
--- /dev/null
+++ b/ITensorGaussianMPS/test/dynamics.jl
@@ -0,0 +1,64 @@
+using ITensorGaussianMPS
+using ITensors
+using LinearAlgebra
+using Test
+
+function diagonalize_h(N)
+  hterms = OpSum()
+  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)
+  ϵ, ϕ = eigen(h)
+  return ϵ, ϕ
+end
+
+@testset "Green Function Properties" begin
+  N = 20
+  ϵ, ϕ = diagonalize_h(N)
+
+  t = 1.0
+
+  Gr = retarded_green_function(t, ϕ, ϵ)
+  Gl = lesser_green_function(t, ϕ, ϵ)
+  Gg = greater_green_function(t, ϕ, ϵ)
+
+  @test Gr ≈ (Gg - Gl)
+
+  # 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
+  ϵ, ϕ = diagonalize_h(N)
+
+  t = 1.0
+
+  G = retarded_green_function(t, ϕ, ϵ)
+
+  G5_10 = retarded_green_function(t, ϕ, ϵ; sites=5:10)
+  @test size(G5_10) == (6, 6)
+  @test G5_10 ≈ G[5:10, 5:10]
+
+  g7 = retarded_green_function(t, ϕ, ϵ; sites=7)
+  @test g7 isa Number
+
+  # Non-contiguous case
+  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
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