Skip to content

Commit

Permalink
[ITensors] [ENHANCMENT] Non-Hermitian dmrg (#913)
Browse files Browse the repository at this point in the history
* nonHermitician

* nonHermitician

* nonHermitician

* nonHermitician

* Change DMRGObserver to allow complex energies

* Change DMRGObserver to allow complex energies

* Formatting and cleanup of tests

* Replace complex_energy with energy_type

* Redesign of DMRGObserver constructors

Co-authored-by: leachinellato <leachinellato@outlook.com.ar>
Co-authored-by: Miles <miles.stoudenmire@gmail.com>
  • Loading branch information
3 people authored May 14, 2022
1 parent 6695289 commit bd31af8
Show file tree
Hide file tree
Showing 4 changed files with 176 additions and 26 deletions.
28 changes: 12 additions & 16 deletions src/mps/dmrg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -29,15 +29,15 @@ and the `sweeps` object determines the parameters used to
control the DMRG algorithm.
Returns:
* `energy::Float64` - eigenvalue of the optimized MPS
* `energy::Complex` - eigenvalue of the optimized MPS
* `psi::MPS` - optimized MPS
Optional keyword arguments:
* `outputlevel::Int = 1` - larger outputlevel values make DMRG print more information and 0 means no output
* `observer` - object implementing the [Observer](@ref observer) interface which can perform measurements and stop DMRG early
* `write_when_maxdim_exceeds::Int` - when the allowed maxdim exceeds this value, begin saving tensors to disk to free memory in large calculations
"""
function dmrg(H::MPO, psi0::MPS, sweeps::Sweeps; kwargs...)::Tuple{Number,MPS}
function dmrg(H::MPO, psi0::MPS, sweeps::Sweeps; kwargs...)
check_hascommoninds(siteinds, H, psi0)
check_hascommoninds(siteinds, H, psi0')
# Permute the indices to have a better memory layout
Expand Down Expand Up @@ -65,10 +65,10 @@ the set of MPOs [H1,H2,H3,..] is efficiently looped over at
each step of the DMRG algorithm when optimizing the MPS.
Returns:
* `energy::Float64` - eigenvalue of the optimized MPS
* `energy::Complex` - eigenvalue of the optimized MPS
* `psi::MPS` - optimized MPS
"""
function dmrg(Hs::Vector{MPO}, psi0::MPS, sweeps::Sweeps; kwargs...)::Tuple{Number,MPS}
function dmrg(Hs::Vector{MPO}, psi0::MPS, sweeps::Sweeps; kwargs...)
for H in Hs
check_hascommoninds(siteinds, H, psi0)
check_hascommoninds(siteinds, H, psi0')
Expand All @@ -95,12 +95,10 @@ and the `sweeps` object determines the parameters used to
control the DMRG algorithm.
Returns:
* `energy::Float64` - eigenvalue of the optimized MPS
* `energy::Complex` - eigenvalue of the optimized MPS
* `psi::MPS` - optimized MPS
"""
function dmrg(
H::MPO, Ms::Vector{MPS}, psi0::MPS, sweeps::Sweeps; kwargs...
)::Tuple{Number,MPS}
function dmrg(H::MPO, Ms::Vector{MPS}, psi0::MPS, sweeps::Sweeps; kwargs...)
check_hascommoninds(siteinds, H, psi0)
check_hascommoninds(siteinds, H, psi0')
for M in Ms
Expand All @@ -113,7 +111,7 @@ function dmrg(
return dmrg(PMM, psi0, sweeps; kwargs...)
end

function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)::Tuple{Number,MPS}
function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)
if length(psi0) == 1
error(
"`dmrg` currently does not support system sizes of 1. You can diagonalize the MPO tensor directly with tools like `LinearAlgebra.eigen`, `KrylovKit.eigsolve`, etc.",
Expand All @@ -137,7 +135,7 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)::Tuple{Number,MPS}
)

# eigsolve kwargs
eigsolve_tol::Float64 = get(kwargs, :eigsolve_tol, 1e-14)
eigsolve_tol::Number = get(kwargs, :eigsolve_tol, 1e-14)
eigsolve_krylovdim::Int = get(kwargs, :eigsolve_krylovdim, 3)
eigsolve_maxiter::Int = get(kwargs, :eigsolve_maxiter, 1)
eigsolve_verbosity::Int = get(kwargs, :eigsolve_verbosity, 0)
Expand Down Expand Up @@ -229,7 +227,8 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)::Tuple{Number,MPS}
maxiter=eigsolve_maxiter,
)
end
energy::Number = vals[1]

energy = vals[1]
phi::ITensor = vecs[1]

ortho = ha == 1 ? "left" : "right"
Expand Down Expand Up @@ -269,9 +268,7 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)::Tuple{Number,MPS}
end

if outputlevel >= 2
@printf(
"Sweep %d, half %d, bond (%d,%d) energy=%.12f\n", sw, ha, b, b + 1, energy
)
@printf("Sweep %d, half %d, bond (%d,%d) energy=%s\n", sw, ha, b, b + 1, energy)
@printf(
" Truncated using cutoff=%.1E maxdim=%d mindim=%d\n",
cutoff(sweeps, sw),
Expand Down Expand Up @@ -300,7 +297,7 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)::Tuple{Number,MPS}
end
if outputlevel >= 1
@printf(
"After sweep %d energy=%.12f maxlinkdim=%d maxerr=%.2E time=%.3f\n",
"After sweep %d energy=%s maxlinkdim=%d maxerr=%.2E time=%.3f\n",
sw,
energy,
maxlinkdim(psi),
Expand All @@ -310,7 +307,6 @@ function dmrg(PH, psi0::MPS, sweeps::Sweeps; kwargs...)::Tuple{Number,MPS}
flush(stdout)
end
isdone = checkdone!(obs; energy=energy, psi=psi, sweep=sw, outputlevel=outputlevel)

isdone && break
end
return (energy, psi)
Expand Down
38 changes: 28 additions & 10 deletions src/mps/observer.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,40 +27,50 @@ implements custom measurements and allows
the `dmrg` function to return early if an
energy convergence criterion is met.
"""
struct DMRGObserver <: AbstractObserver
struct DMRGObserver{T} <: AbstractObserver
ops::Vector{String}
sites::Vector{<:Index}
measurements::Dict{String,DMRGMeasurement}
energies::Vector{Float64}
energies::Vector{T}
truncerrs::Vector{Float64}
etol::Float64
minsweeps::Int64
end

"""
DMRGObserver(;energy_tol=0.0,
minsweeps=2)
minsweeps=2,
energy_type=Float64)
Construct a DMRGObserver by providing the energy
tolerance used for early stopping, and minimum number
of sweeps that must be done.
Optional keyword arguments:
- energy_tol: if the energy from one sweep to the
next no longer changes by more than this amount,
stop after the current sweep
- minsweeps: do at least this many sweeps
- energy_type: type to use when storing energies at each step
"""
function DMRGObserver(; energy_tol=0.0, minsweeps=2)
function DMRGObserver(; energy_tol=0.0, minsweeps=2, energy_type=Float64)
return DMRGObserver(
[], Index[], Dict{String,DMRGMeasurement}(), [], [], energy_tol, minsweeps
String[],
Index[],
Dict{String,DMRGMeasurement}(),
energy_type[],
Float64[],
energy_tol,
minsweeps,
)
end

"""
DMRGObserver(ops::Vector{String},
sites::Vector{<:Index};
energy_tol=0.0,
minsweeps=2)
minsweeps=2,
energy_type=Float64)
Construct a DMRGObserver, provide an array
of `ops` of operator names which are strings
Expand All @@ -76,16 +86,24 @@ Optionally, one can provide an energy
tolerance used for early stopping, and minimum number
of sweeps that must be done.
Optional keyword arguments:
- energy_tol: if the energy from one sweep to the
next no longer changes by more than this amount,
stop after the current sweep
- minsweeps: do at least this many sweeps
- energy_type: type to use when storing energies at each step
"""
function DMRGObserver(
ops::Vector{String}, sites::Vector{<:Index}; energy_tol=0.0, minsweeps=2
ops::Vector{String},
sites::Vector{<:Index};
energy_tol=0.0,
minsweeps=2,
energy_type=Float64,
)
measurements = Dict(o => DMRGMeasurement() for o in ops)
return DMRGObserver(ops, sites, measurements, [], [], energy_tol, minsweeps)
return DMRGObserver{energy_type}(
ops, sites, measurements, energy_type[], Float64[], energy_tol, minsweeps
)
end

"""
Expand Down Expand Up @@ -157,8 +175,8 @@ end
function checkdone!(o::DMRGObserver; kwargs...)
outputlevel = get(kwargs, :outputlevel, false)
if (
length(energies(o)) > o.minsweeps &&
abs(energies(o)[end] - energies(o)[end - 1]) < o.etol
length(real(energies(o))) > o.minsweeps &&
abs(real(energies(o))[end] - real(energies(o))[end - 1]) < o.etol
)
outputlevel > 0 && println("Energy difference less than $(o.etol), stopping DMRG")
return true
Expand Down
68 changes: 68 additions & 0 deletions test/XXZ_complex.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
#########
## XXZ model with imaginary part
#########
using ITensors

function Hamiltonian(sites, Δ::Float64, J::Float64, γ::Float64, h::Float64)
N = length(sites)
ampo = AutoMPO()

for j in 1:(N - 1)
ampo += Δ, "Sz", j, "Sz", j + 1
ampo += -(J + im * γ * (-1.0)^j), "Sx", j, "Sx", j + 1
ampo += -(J + im * γ * (-1.0)^j), "Sy", j, "Sy", j + 1
end
for j in 1:N
ampo += h * (-1.0)^j, "Sz", j
end
# Convert these terms to an MPO tensor network
return MPO(ampo, sites)
end

let
#model parameters
N = 4
Δ = -0.70
J = 1.0
γ = 0.1
h = 0.1

alg = "qr_iteration"

#dmrg parameters
sweeps = Sweeps(1000)
minsweeps = 5
maxdim!(sweeps, 50, 100, 200)
#cutoff!(sweeps, 1E-12)
etol = 1E-12

sites = siteinds("S=1/2", N; conserve_qns=false)

#initial state
state = ["Emp" for n in 1:N]
p = N
for i in N:-1:1
if p > i
#println("Doubly occupying site $i")
state[i] = "UpDn"
p -= 2
elseif p > 0
#println("Singly occupying site $i")
state[i] = (isodd(i) ? "Up" : "Dn")
p -= 1
end
end
psi0 = randomMPS(sites, state)
@show flux(psi0)

H = Hamiltonian(sites, Δ, J, γ, h)

obs = DMRGObserver(
["Sz"], sites; energy_tol=etol, minsweeps=minsweeps, complex_energies=true
)

energy, psi = dmrg(
H, psi0, sweeps; svd_alg=alg, observer=obs, outputlevel=1, ishermitian=false
)
println("Final energy = $energy")
end
68 changes: 68 additions & 0 deletions test/trans_ising_complex.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@
using ITensors

function Hamiltonian(sites, λ::Float64, k::Float64)
N = length(sites)
ampo = AutoMPO()

for j in 1:(N - 1)
ampo += -0.5 * 2 * λ, "Sz", j, "Sz", j + 1
end
for j in 1:N
ampo += -0.5 * 2 * im * k, "Sz", j
ampo += -0.5, "S+", j
ampo += -0.5, "S-", j
end
# Convert these terms to an MPO tensor network
return MPO(ampo, sites)
end

let
#-----------------------------------------------------------------------

#models parameters
N = 4
λ = 0.1
k = 0.5

alg = "qr_iteration"

#dmrg parameters
sweeps = Sweeps(1000)
minsweeps = 5
maxdim!(sweeps, 50, 100, 200)
#cutoff!(sweeps, 1E-12)
etol = 1E-12

#-----------------------------------------------------------------------

sites = siteinds("S=1/2", N; conserve_qns=false)
#-----------------------------------------------------------------------

#inicial state
state = ["Emp" for n in 1:N]
p = N
for i in N:-1:1
if p > i
state[i] = "UpDn"
p -= 2
elseif p > 0
state[i] = (isodd(i) ? "Up" : "Dn")
p -= 1
end
end
psi0 = randomMPS(sites, state)
@show flux(psi0)
#-----------------------------------------------------------------------

#-----------------------------------------------------------------------

H = Hamiltonian(sites, λ, k)
#-----------------------------------------------------------------------

obs = DMRGObserver(; energy_tol=etol, minsweeps=2, complex_energies=true)

energy, psi = dmrg(
H, psi0, sweeps; svd_alg=alg, observer=obs, outputlevel=1, ishermitian=false
)
println("Final energy = $energy")
end

0 comments on commit bd31af8

Please sign in to comment.