Algorithms for MPS#
DMRG#
In order to execute DMRG we need the hermitian operator (typically a Hamiltonian) written as MPO and an initial guess for MPS. Here is a simple example of DMRG used to obtain the ground state of quadratic Hamiltonian:
def dmrg_XX_model_dense(config=None, tol=1e-6):
"""
Initialize random MPS of dense tensors and
runs a few sweeps of DMRG with the Hamiltonian of XX model.
"""
# Knowing the exact solution we can compare it with the DMRG result.
# In this test we will consider sectors of different occupations.
#
# In this example we use yastn.Tensor with no symmetry imposed.
#
# Here, the Hamiltonian for N = 7 sites
# is obtained with the automatic generator.
#
N = 7
#
opts_config = {} if config is None else \
{'backend': config.backend,
'default_device': config.default_device}
# pytest uses config to inject various backends and devices for testing
ops = yastn.operators.Spin12(sym='dense', **opts_config)
generate = mps.Generator(N=N, operators=ops)
parameters = {"t": 1.0, "mu": 0.2,
"rN": range(N),
"rNN": [(i, i+1) for i in range(N - 1)]}
H_str = r"\sum_{i,j \in rNN} t ( sp_{i} sm_{j} + sp_{j} sm_{i} )"
H_str += r" + \sum_{j \in rN} mu sp_{j} sm_{j}"
H = generate.mpo_from_latex(H_str, parameters)
#
# and MPO to measure occupation:
#
O_occ = generate.mpo_from_latex(r"\sum_{j \in rN} sp_{j} sm_{j}",
parameters)
#
# Known energies and occupations of low-energy eigenstates.
#
Eng_gs = [-3.427339492125, -3.227339492125, -2.861972627395]
occ_gs = [3, 4, 2]
#
# To standardize this test we fix a seed for random MPS we use.
#
generate.random_seed(seed=0)
#
# Set options for truncation for '2site' method of mps.dmrg_.
#
Dmax = 8
opts_svd = {'tol': 1e-8, 'D_total': Dmax}
#
# Finally run DMRG starting from random MPS psi
#
psi = generate.random_mps(D_total=Dmax)
#
# A single run for a ground state can be done using:
# mps.dmrg_(psi, H, method=method,energy_tol=tol/10,
# max_sweeps=20, opts_svd=opts_svd)
#
# We create a subfunction run_dmrg to explain how to
# target some sectors for occupation. This is not necessary
# but we do it for the sake of clarity and testing.
#
psis = run_dmrg(psi, H, O_occ, Eng_gs, occ_gs, opts_svd, tol)
#
# _dmrg can be executed as a generator to monitor states
# between dmrg sweeps. This is done by providing `iterator_step`.
#
psi = generate.random_mps(D_total=Dmax)
for step in mps.dmrg_(psi, H, method='1site',
max_sweeps=20, iterator_step=2):
assert step.sweeps % 2 == 0 # stop every second sweep here
occ = mps.measure_mpo(psi, O_occ, psi) # measure occupation
# here breaks if it is close to the known result.
if abs(occ.item() - occ_gs[0]) < tol:
break
assert step.sweeps < 20
def run_dmrg(phi, H, O_occ, E_target, occ_target, opts_svd, tol):
r"""
Run mps.dmrg_ to find the ground state and
a few low-energy states of the Hamiltonian H.
Verify resulting energies against known reference solutions.
"""
#
# DMRG can look for solutions in a space
# orthogonal to provided MPSs. We start with empty list,
# project = [], and keep adding to it previously found eigenstates.
# This allows us to find the ground state
# and a few consequative lowest-energy eigenstates.
#
project, states = [], []
for ref_eng, ref_occ in zip(E_target, occ_target):
#
# We find a state and check that its energy and total occupation
# matches the expected reference values.
#
# We copy initial random MPS psi, as
# yastn.dmrg_ modifies provided input state in place.
#
psi = phi.shallow_copy()
#
# We set up dmrg_ to terminate iterations
# when energy is converged within some tolerance.
#
out = mps.dmrg_(psi, H, project=project, method='2site',
energy_tol=tol / 10, max_sweeps=20, opts_svd=opts_svd)
#
# Output of _dmrg is a nametuple with information about the run,
# including the final energy.
# Occupation number has to be calcuted using measure_mpo.
#
eng = out.energy
occ = mps.measure_mpo(psi, O_occ, psi)
#
# Print the result:
#
print(f"2site DMRG; energy: {eng:{1}.{8}} / {ref_eng:{1}.{8}}; "
+ f"occupation: {occ:{1}.{8}} / {ref_occ}")
assert abs(eng - ref_eng) < tol * 100
assert abs(occ - ref_occ) < tol
#
# We further iterate with '1site' DMRG
# and stricter convergence criterion.
#
out = mps.dmrg_(psi, H, project=project, method='1site',
Schmidt_tol=tol / 10, max_sweeps=20)
eng = mps.measure_mpo(psi, H, psi)
occ = mps.measure_mpo(psi, O_occ, psi)
print(f"1site DMRG; energy: {eng:{1}.{8}} / {ref_eng:{1}.{8}}; "
+ f"occupation: {occ:{1}.{8}} / {ref_occ}")
# test that energy outputed by dmrg is correct
assert abs(eng - ref_eng) < tol
assert abs(occ - ref_occ) < tol
assert abs(eng - out.energy) < tol # test dmrg_ output information
#
# Finally, we add the found state psi to the list of states
# to be projected out in the next step of the loop.
#
penalty = 100
project.append((penalty, psi))
states.append(psi)
return states
The same can be done for other symmetries:
def dmrg_XX_model_Z2(config=None, tol=1e-6):
"""
Initialize random MPS of Z2 tensors and tests mps.dmrg_ vs known results.
"""
opts_config = {} if config is None else \
{'backend': config.backend,
'default_device': config.default_device}
# pytest uses config to inject various backends and devices for testing
ops = yastn.operators.SpinlessFermions(sym='Z2', **opts_config)
generate = mps.Generator(N=7, operators=ops)
generate.random_seed(seed=0)
N, Dmax = 7, 8
opts_svd = {'tol': 1e-8, 'D_total': Dmax}
Eng_occ_target = {
0: ([-3.227339492125, -2.861972627395, -2.461972627395],
[4, 2, 4]),
1: ([-3.427339492125, -2.661972627395, -2.261972627395],
[3, 3, 5])}
H_str = r"\sum_{i,j \in rNN} t (cp_{i} c_{j} + cp_{j} c_{i})"
H_str += r" + \sum_{j\in rN} mu cp_{j} c_{j}"
parameters = {"t": 1.0, "mu": 0.2,
"rN": range(N),
"rNN": [(i, i+1) for i in range(N - 1)]}
H = generate.mpo_from_latex(H_str, parameters)
O_occ = generate.mpo_from_latex(r"\sum_{j\in rN} cp_{j} c_{j}", parameters)
for parity, (E_target, occ_target) in Eng_occ_target.items():
psi = generate.random_mps(D_total=Dmax, n=parity)
run_dmrg(psi, H, O_occ, E_target, occ_target, opts_svd, tol)
def dmrg_XX_model_U1(config=None, tol=1e-6):
"""
Initialize random MPS of U1 tensors and tests _dmrg vs known results.
"""
opts_config = {} if config is None else \
{'backend': config.backend,
'default_device': config.default_device}
# pytest uses config to inject various backends and devices for testing
ops = yastn.operators.SpinlessFermions(sym='U1', **opts_config)
generate = mps.Generator(N=7, operators=ops)
generate.random_seed(seed=0)
N = 7
H_str = r"\sum_{i,j \in rNN} t (cp_{i} c_{j} + cp_{j} c_{i})"
H_str += r" + \sum_{j\in rN} mu cp_{j} c_{j}"
parameters = {"t": 1.0, "mu": 0.2,
"rN": range(N),
"rNN": [(i, i+1) for i in range(N - 1)]}
H = generate.mpo_from_latex(H_str, parameters)
O_occ = generate.mpo_from_latex(r"\sum_{j\in rN} cp_{j} c_{j}", parameters)
Eng_sectors = {
2: [-2.861972627395, -2.213125929752], # -1.779580427103],
3: [-3.427339492125, -2.661972627395], # -2.013125929752],
4: [-3.227339492125, -2.461972627395], # -1.813125929752]
}
Dmax = 8
opts_svd = {'tol': 1e-8, 'D_total': Dmax}
for occ_sector, E_target in Eng_sectors.items():
psi = generate.random_mps(D_total=Dmax, n=occ_sector)
occ_target = [occ_sector] * len(E_target)
run_dmrg(psi, H, O_occ, E_target, occ_target, opts_svd, tol)
See as well the examples for Multiplication, which contains DMRG and variational MPS compression.
TDVP#
Sudden quench in a free-fermionic model#
def tdvp_sudden_quench(sym='U1', config=None, tol=1e-10):
"""
Simulate a sudden quench of a free-fermionic (hopping) model.
Compare observables versus known reference results.
"""
N, n = 6, 3 # Consider a system of 6 modes and 3 particles.
#
# Load operators
#
opts_config = {} if config is None else \
{'backend': config.backend,
'default_device': config.default_device}
# pytest uses config to inject various backends and devices for testing
ops = yastn.operators.SpinlessFermions(sym=sym, **opts_config)
ops.random_seed(seed=0)
#
# Hopping matrix, functions consuming it use only upper-triangular part.
#
J0 = [[1, 0.5j, 0, 0.3, 0.1, 0 ],
[0, -1, 0.5j, 0, 0.3, 0.1 ],
[0, 0, 1, 0.5j, 0, 0.3 ],
[0, 0, 0, -1, 0.5j, 0 ],
[0, 0, 0, 0, 1, 0.5j],
[0, 0, 0, 0, 0, -1 ]]
#
# Generate corresponding MPO using the function from previous examples
#
H0 = build_mpo_hopping_Hterm(J0, sym=sym, config=config)
#
# Find the ground state using DMRG
# Bond dimension Dmax = 8 for N = 6 is large enough
# to avoid truncation errors in the test
#
Dmax = 8
opts_svd = {'tol': 1e-15, 'D_total': Dmax}
#
# Run DMRG for the ground state
# global ground state is for n=3, i.e., works for Z2 and U1
#
I = mps.product_mpo(ops.I(), N)
n_psi = n % 2 if sym=='Z2' else n # for U1; charge of MPS
psi = mps.random_mps(I, D_total=Dmax, n=n_psi)
#
out = mps.dmrg_(psi, H0, method='2site', max_sweeps=2, opts_svd=opts_svd)
out = mps.dmrg_(psi, H0, method='1site', max_sweeps=10,
energy_tol=1e-14, Schmidt_tol=1e-14)
#
# Get reference results for the ground state and check mps
#
C0ref, E0ref = gs_correlation_matrix_exact(J0, n) # defined below
assert abs(out.energy - E0ref) < tol
#
# Calculate correlation matrix for MPS and test vs exact reference
#
C0psi = correlation_matrix_from_mps(psi, ops, tol) # defined below
assert np.allclose(C0ref, C0psi, rtol=tol)
#
# Sudden quench with a new Hamiltonian
#
J1 = [[-1, 0.5, 0, -0.3, 0.1, 0 ],
[ 0, 1 , 0.5, 0, -0.3, 0.1],
[ 0, 0 , -1, 0.5, 0, -0.3],
[ 0, 0 , 0, 1, 0.5, 0 ],
[ 0, 0 , 0, 0, -1, 0.5],
[ 0, 0 , 0, 0, 0, 1 ]]
H1 = build_mpo_hopping_Hterm(J1, sym=sym, config=config)
#
# Run time evolution and calculate correlation matrix at two snapshots
#
times = (0, 0.25, 0.6)
#
# Parameters for expmv in tdvp_,
# 'ncv' is an initial guess for the size of Krylov space.
# It gets updated at each site/bond during evolution.
#
opts_expmv = {'hermitian': True, 'ncv': 5, 'tol': 1e-12}
#
for method in ('1site', '2site', '12site'): # test various methods
# shallow_copy is sufficient to retain the initial state
phi = psi.shallow_copy()
for step in mps.tdvp_(phi, H1, times=times, method=method, dt=0.125,
opts_svd=opts_svd, opts_expmv=opts_expmv):
#
# Calculate correlation matrix from mps.
# Compare with exact reference results defined below.
#
Cphi = correlation_matrix_from_mps(phi, ops, tol)
Cref = evolve_correlation_matrix_exact(C0ref, J1, step.tf)
assert np.allclose(Cref, Cphi, rtol=tol)
def correlation_matrix_from_mps(psi, ops, tol):
"""
Calculate correlation matrix for MPS psi.
"""
assert abs(psi.norm() - 1) < tol # check normalization
cpc = mps.measure_2site(psi, ops.cp(), ops.c(), psi, bonds='<=>') # all
C = np.zeros((psi.N, psi.N), dtype=np.complex128)
for (n1, n2), v in cpc.items():
C[n2, n1] = v
return C
def gs_correlation_matrix_exact(J, n):
"""
Correlation matrix for the ground state
of n particles with hopping Hamiltonian matrix J.
C[m, n] = <c_n^dag c_m>
"""
J = np.triu(J, 0) + np.triu(J, 1).T.conj()
D, V = np.linalg.eigh(J)
Egs = np.sum(D[:n])
C0 = np.zeros(len(D))
C0[:n] = 1
C = V @ np.diag(C0) @ V.T.conj()
return C, Egs
def evolve_correlation_matrix_exact(Ci, J, t):
"""
Evolve correlation matrix C by time t with hopping Hamiltonian J.
Diagonal elements of J array is on-site potential and
upper triangular terms of J are hopping amplitudes.
"""
J = np.triu(J, 0) + np.triu(J, 1).T.conj()
# U = expm(1j * t * J)
D, V = np.linalg.eigh(J)
U = V @ np.diag(np.exp(1j * t * D)) @ V.T.conj()
Cf = U.conj().T @ Ci @ U
return Cf
Slow quench across a quantum critical point in a transverse Ising chain#
def tdvp_KZ_quench(sym='Z2', config=None):
"""
Simulate a slow quench across a quantum critical point in
a small transverse field Ising chain with periodic boundary conditions.
Compare with exact reference results.
"""
#
N = 10 # Consider a system of 10 sites
#
# Load spin-1/2 operators
#
opts_config = {} if config is None else \
{'backend': config.backend,
'default_device': config.default_device}
# pytest uses config to inject various backends and devices for testing
#
ops = yastn.operators.Spin12(sym=sym, **opts_config)
ops.random_seed(seed=0)
#
# Hterm-s to generate H = -sum_i X_i X_{i+1} - g * Z_i
#
I = mps.product_mpo(ops.I(), N) # identity MPO
termsXX = [mps.Hterm(-1, [i, (i + 1) % N], [ops.x(), ops.x()]) for i in range(N)]
HXX = mps.generate_mpo(I, termsXX)
termsZ = [mps.Hterm(-1, [i], [ops.z()]) for i in range(N)]
HZ = mps.generate_mpo(I, termsZ)
#
# Kibble-Zurek quench across a critical point at gc = 1
# tauQ is the quench time
#
tauQ, gc = 1, 1
ti, tf = -tauQ, tauQ # evolve from gi = 2 to gf = 0
H = lambda t: [HXX, (gc - t / tauQ) * HZ] # linear quench
#
# Analytical reference expectation values measured
# at g = 1 and g = 0 (for tauQ=1, gi=2 and N=10)
#
XXex = {1: 0.470182292934, 0: 0.738769410121}
Zex = {1: 0.776255260472, 0: 0.163491011822}
Egs = -2.127120881869 # the ground state energy at gi = 2
#
# Start with the ground state at gi = 2
# Run DMRG to get the initial ground state
#
Dmax = 10
psi = mps.random_mps(I, D_total=Dmax)
out = mps.dmrg_(psi, H(ti), method='2site', max_sweeps=2,
opts_svd={'tol': 1e-6, 'D_total': Dmax})
out = mps.dmrg_(psi, H(ti), method='1site', max_sweeps=10,
energy_tol=1e-12, Schmidt_tol=1e-12)
#
# Test ground state energy versus reference
#
assert abs(out.energy / N - Egs) < 1e-7
assert psi.get_bond_dimensions() == (1, 2, 4, 8, 10, 10, 10, 8, 4, 2, 1)
#
# Slow quench to gf = 0
# Sets up tdvp_ parameters; allows the growth of the bond dimension.
#
Dmax = 16
opts_expmv = {'hermitian': True, 'tol': 1e-12}
opts_svd = {'tol': 1e-6, 'D_total': Dmax}
#
for step in mps.tdvp_(psi, H, times=(ti, 0, tf),
method='12site',dt=0.04, order='2nd',
opts_svd=opts_svd, opts_expmv=opts_expmv):
#
# tdvp_() always gives an iterator
# Calculate expectation values at snapshots
#
EZ = mps.measure_1site(psi, ops.z(), psi)
EXX = mps.measure_2site(psi, ops.x(), ops.x(), psi, bonds='r1p') # periodic nn
#
# Compare them with the exact result
#
gg = round(gc - step.tf / tauQ) # g at the snapshot
assert all(abs(EXX[k, (k + 1) % N] - XXex[gg]) < 1e-4 for k in range(N))
assert all(abs(EZ[k] - Zex[gg]) < 1e-4 for k in range(N))
#
# Bond dimension was updated
#
bd_ref = (1, 2, 4, 8, 16, 16, 16, 8, 4, 2, 1)
assert psi.get_bond_dimensions() == bd_ref