Initialization and creation operations#
Building MPS/MPO manually#
The tensor making up MPS/MPO can be assigned manually, setting them one by one.
Note that the virtual dimensions/spaces of the neighboring MPS/MPO tensors should be consistent.
Ground state of spin-1 AKLT model#
Here, as an example, we set up a well-known exact MPS: A ground state of Affleck-Kennedy-Lieb-Tasaki (AKLT) model.
import numpy as np
import yastn
import yastn.tn.mps as mps
def build_aklt_state_manually(N=5, lvec=(1, 0), rvec=(0, 1), config=None):
"""
Initialize MPS tensors by hand. Example for Spin-1 AKLT state of N sites.
Allow inputing boundary vectors for an MPS with open boundary conditions.
"""
# Prepare rank-3 on-site tensor with virtual dimensions 2
# and physical dimension dim(Spin-1)=3
# _
# dim(left)=2 --|A|-- dim(right)=2
# |
# dim(Spin-1)=3
#
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
config = yastn.make_config(**opts_config)
A = yastn.Tensor(config, s=(-1, 1, 1))
s13, s23 = np.sqrt(1. / 3), np.sqrt(2. / 3)
A.set_block(Ds=(2, 3, 2), val=[[[0, s23], [-s13, 0], [0, 0]],
[[0, 0], [0, s13], [-s23, 0]]])
#
# as well as left and right boundary vectors
#
LV = yastn.Tensor(config, s=(-1, 1))
LV.set_block(Ds=(1, 2), val=lvec)
RV = yastn.Tensor(config, s=(-1, 1))
RV.set_block(Ds=(2, 1), val=rvec)
#
# We initialize empty MPS for N sites
# and assign its on-site tensors one-by-one.
#
psi = mps.Mps(N)
for n in range(N):
psi[n] = A.copy()
#
# Due to open boundary conditions, the first and the last
# MPS tensors have left and right virtual indices projected to dimension 1
#
psi[psi.first] = LV @ psi[psi.first]
psi[psi.last] = psi[psi.last] @ RV
#
# The resulting MPS is not normalized, nor in any canonical form.
#
return psi
Hamiltonian for nearest-neighbor hopping/XX model#
The same can be done for MPOs. Here, we show a construction of a simple nearest-neighbour hopping Hamiltonian with hopping amplitude \(t\) and on-site energy \(\mu\) with different realizations of explicit symmetry.
def build_mpo_nn_hopping_manually(N, t, mu, sym='U1', config=None):
"""
Nearest-neighbor hopping Hamiltonian on N sites
with hopping amplitude t and chemical potential mu,
H = t * sum_{n=1}^{N-1} (cp_n c_{n+1} + cp_{n+1} c_n)
+ mu * sum_{n=1}^{N} cp_n c_n
Initialize MPO symmetric tensors by hand with sym = 'dense', 'Z2', or 'U1'.
Config is used to inject non-default backend and default_device.
"""
#
# Build empty MPO for system of N sites
#
H = mps.Mpo(N)
#
# Depending on the symmetry, define elements of on-site tensor
#
# We chose signature convention for indices of the MPO tensor as follows
# |
# ^(-1)
# |
# (-1)-<-|T|-<-(+1)
# |
# ^(+1)
# |
#
# We set fermionic=True for consistency as in tests
# MPO is compared with operators.SpinlessFermions.
# backend and device is inherited from config to automitize pytest testing.
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
config = yastn.make_config(fermionic=True, sym=sym, **opts_config)
if sym == 'dense': # no symmetry
# Basic rank-2 blocks (matrices) of on-site tensors
cp = np.array([[0, 0], [1, 0]])
c = np.array([[0, 1], [0, 0]])
nn = np.array([[0, 0], [0, 1]])
ee = np.array([[1, 0], [0, 1]])
oo = np.array([[0, 0], [0, 0]])
for n in H.sweep(to='last'):
# empty tensors
H[n] = yastn.Tensor(config=config, s=(-1, 1, 1, -1))
# that gets filled-in
if n == H.first:
tmp = np.block([[mu * nn, t * cp, t * c, ee]])
H[n].set_block(val=tmp, Ds=(1, 2, 4, 2))
elif n == H.last:
tmp = np.block([[ee], [c], [cp], [mu * nn]])
H[n].set_block(val=tmp, Ds=(4, 2, 1, 2))
else:
tmp = np.block([[ee, oo, oo, oo],
[c, oo, oo, oo],
[cp, oo, oo, oo],
[mu * nn, t * cp, t * c, ee]])
H[n].set_block(val=tmp, Ds=(4, 2, 4, 2))
elif sym == 'Z2': # Z2 symmetry
for n in H.sweep(to='last'):
H[n] = yastn.Tensor(config=config, s=(-1, 1, 1, -1), n=0)
if n == H.first:
H[n].set_block(ts=(0, 0, 0, 0),
Ds=(1, 1, 2, 1), val=[0, 1])
H[n].set_block(ts=(0, 1, 0, 1),
Ds=(1, 1, 2, 1), val=[mu, 1])
H[n].set_block(ts=(0, 0, 1, 1),
Ds=(1, 1, 2, 1), val=[t, 0])
H[n].set_block(ts=(0, 1, 1, 0),
Ds=(1, 1, 2, 1), val=[0, t])
elif n == H.last:
H[n].set_block(ts=(0, 0, 0, 0),
Ds=(2, 1, 1, 1), val=[1, 0])
H[n].set_block(ts=(0, 1, 0, 1),
Ds=(2, 1, 1, 1), val=[1, mu])
H[n].set_block(ts=(1, 1, 0, 0),
Ds=(2, 1, 1, 1), val=[1, 0])
H[n].set_block(ts=(1, 0, 0, 1),
Ds=(2, 1, 1, 1), val=[0, 1])
else:
H[n].set_block(ts=(0, 0, 0, 0),
Ds=(2, 1, 2, 1), val=[[1, 0], [0, 1]])
H[n].set_block(ts=(0, 1, 0, 1),
Ds=(2, 1, 2, 1), val=[[1, 0], [mu, 1]])
H[n].set_block(ts=(0, 0, 1, 1),
Ds=(2, 1, 2, 1), val=[[0, 0], [t, 0]])
H[n].set_block(ts=(0, 1, 1, 0),
Ds=(2, 1, 2, 1), val=[[0, 0], [0, t]])
H[n].set_block(ts=(1, 1, 0, 0),
Ds=(2, 1, 2, 1), val=[[1, 0], [0, 0]])
H[n].set_block(ts=(1, 0, 0, 1),
Ds=(2, 1, 2, 1), val=[[0, 0], [1, 0]])
elif sym == 'U1': # U1 symmetry
for n in H.sweep(to='last'):
H[n] = yastn.Tensor(config=config, s=(-1, 1, 1, -1), n=0)
if n == H.first:
H[n].set_block(ts=(0, 0, 0, 0),
val=[0, 1], Ds=(1, 1, 2, 1))
H[n].set_block(ts=(0, 1, 0, 1),
val=[mu, 1], Ds=(1, 1, 2, 1))
H[n].set_block(ts=(0, 0, 1, 1),
val=[t], Ds=(1, 1, 1, 1))
H[n].set_block(ts=(0, 1, -1, 0),
val=[t], Ds=(1, 1, 1, 1))
elif n == H.last:
H[n].set_block(ts=(0, 0, 0, 0),
Ds=(2, 1, 1, 1), val=[1, 0])
H[n].set_block(ts=(0, 1, 0, 1),
Ds=(2, 1, 1, 1), val=[1, mu])
H[n].set_block(ts=(1, 1, 0, 0),
Ds=(1, 1, 1, 1), val=[1])
H[n].set_block(ts=(-1, 0, 0, 1),
Ds=(1, 1, 1, 1), val=[1])
else:
H[n].set_block(ts=(0, 0, 0, 0),
Ds=(2, 1, 2, 1), val=[[1, 0], [0, 1]])
H[n].set_block(ts=(0, 1, 0, 1),
Ds=(2, 1, 2, 1), val=[[1, 0], [mu, 1]])
H[n].set_block(ts=(0, 0, 1, 1),
Ds=(2, 1, 1, 1), val=[0, t])
H[n].set_block(ts=(0, 1, -1, 0),
Ds=(2, 1, 1, 1), val=[0, t])
H[n].set_block(ts=(1, 1, 0, 0),
Ds=(1, 1, 2, 1), val=[1, 0])
H[n].set_block(ts=(-1, 0, 0, 1),
Ds=(1, 1, 2, 1), val=[1, 0])
return H
Building MPO using Hterm#
Hamiltonian is typically given as a sum of products of local operators.
We can specify a single such product using mps.tn.mps.Hterm.
Hamiltonian is generated by a function yastn.tn.mps.generate_mpo which takes an identity MPO and a Sequence[Hterm].
Spinless fermions with long-range hopping (Hterm)#
def build_mpo_hopping_Hterm(J, sym="U1", config=None):
"""
Fermionic hopping Hamiltonian on N sites with hoppings at arbitrary range.
The upper triangular part of N x N matrix J defines hopping amplitudes,
and the diagonal defines on-site chemical potentials.
"""
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)
#
Hterms = [] # list of Hterm(amplitude, positions, operators)
#
# Each Hterm corresponds to a single product of local operators.
# Hamiltonian is a sum of such products.
#
N = len(J)
#
# chemical potential on site n
#
for n in range(N):
if abs(J[n][n]) > 0:
Hterms.append(mps.Hterm(amplitude=J[n][n],
positions=[n],
operators=[ops.n()]))
#
# hopping term between sites m and n
for m in range(N):
for n in range(m + 1, N):
if abs(J[m][n]) > 0:
Hterms.append(mps.Hterm(amplitude=J[m][n],
positions=(m, n),
operators=(ops.cp(), ops.c())))
Hterms.append(mps.Hterm(amplitude=np.conj(J[m][n]),
positions=(n, m),
operators=(ops.cp(), ops.c())))
#
# We need an identity MPO operator.
#
I = mps.product_mpo(ops.I(), N)
#
# Generate MPO for Hterms
#
H = mps.generate_mpo(I, Hterms)
#
return H
Generator class for MPO/MPS#
The MPO can be constructed automatically using dedicated yastn.tn.mps.Generator supplied with the set of local operators.
Automatic generator creates MPO with symmetries inherited from local operators.
Generator can be also used to create random MPS and MPO, where local Hilbert spaces are read from identity operators defined in the Generator.
Building random MPS/MPO#
def random_mps_spinless_fermions(N=10, D_total=16, sym='Z2', n=1, config=None):
"""
Generate random MPS of N sites, with bond dimension D_total,
tensors with symmetry sym and total charge n.
"""
#
# load local 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) # fix seed
I = mps.product_mpo(ops.I(), N) # identity MPS
#
# random MPS of matching physical dimension
#
psi = mps.random_mps(I, D_total=D_total, n=n)
#
# the same can be done employing Generator class
#
generate = mps.Generator(N, ops)
generate.random_seed(seed=0)
phi = generate.random_mps(D_total=D_total, n=n)
assert (psi - phi).norm() < 1e-12 * psi.norm()
#
return psi
Hamiltonian for nearest-neighbor hopping/XX model (Generator)#
def mpo_nn_hopping_latex(N, t, mu, sym="U1", config=None):
"""
Nearest-neighbor hopping Hamiltonian on N sites
with hopping amplitude t and chemical potential mu.
"""
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)
Hstr = r"\sum_{j,k \in NN} t (cp_{j} c_{k}+cp_{k} c_{j})"
Hstr += r" + \sum_{i \in sites} mu cp_{i} c_{i}"
parameters = {"t": t,
"mu": mu,
"sites": list(range(N)),
"NN": list((i, i+1) for i in range(N-1))}
generate = mps.Generator(N, ops)
H = generate.mpo_from_latex(Hstr, parameters=parameters)
return H
Save and load MPS/MPO#
MPS/MPO can be saved/loaded either to/from a dictionary or an HDF5 file. YASTN configuration with matching symmetry has to be provided for on-site MPS/MPO tensors when loading either from a dictionary or HDF5 file.
Using Python’s dictionary#
def save_load_mps_dict(sym='dense', config=None, tol=1e-12):
"""
Initialize random MPS and checks saving/loading to/from npy file.
"""
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
#
# generate random mps with 3-dimensional local spaces
#
ops = yastn.operators.Spin1(sym=sym, **opts_config)
I = mps.product_mpo(ops.I(), N=31)
psi = -0.5 * mps.random_mps(I, D_total=25, dtype='complex128')
#
# Next, we serialize MPS into dictionary.
#
tmp = psi.save_to_dict()
#
# Last, we load the MPS from the dictionary,
# providing valid YASTN configuration
#
config = ops.config
phi = mps.load_from_dict(config, tmp)
#
# Test psi == phi
#
assert (psi - phi).norm() < tol * psi.norm()
#
# Similarly for MPO
#
psi = -1j * mps.random_mpo(I, D_total=11) # adding extra complex factor
psi.canonize_(to='last', normalize=False) # extra cannonization
psi.canonize_(to='first', normalize=False) # retaining the norm
tmp = psi.save_to_dict()
phi = mps.load_from_dict(config, tmp)
assert (psi - phi).norm() < tol * psi.norm()
Using HDF5 format#
import os
import h5py
def save_load_mps_hdf5(sym='dense', config=None, tol=1e-12):
"""
Initialize random MPS and checks saving/loading to/from HDF5 file.
"""
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
#
# generate random mps with 3-dimensional local spaces
#
ops = yastn.operators.Spin1(sym=sym, **opts_config)
I = mps.product_mpo(ops.I(), N=31)
psi = 2 * mps.random_mps(I, D_total=25) # adding extra factor
#
# We delete file if it already exists.
# (It is enough to clear the address './state' if the file already exists)
#
try:
os.remove("tmp.h5")
except OSError:
pass
#
# We save the MPS to file 'tmp.h5' under the address 'state/'
#
with h5py.File('tmp.h5', 'a') as f:
psi.save_to_hdf5(f, 'state/')
#
# To read MPS from HDF5 file, open the file and load the MPS stored
# at address 'state/'.
# Note: You have to provide valid YASTN configuration
#
config = ops.config
with h5py.File('tmp.h5', 'r') as f:
phi = mps.load_from_hdf5(config, f, './state/')
#
# Test psi == phi
#
assert (psi - phi).norm() < tol * psi.norm()
#
# Similarily, one can save and load MPO
#
psi = -1j * mps.random_mpo(I, D_total=8, dtype='complex128')
psi.canonize_(to='last', normalize=False) # extra cannonization
psi.canonize_(to='first', normalize=False) # retaining the norm
with h5py.File('tmp.h5', 'w') as f:
psi.save_to_hdf5(f, 'state/')
with h5py.File('tmp.h5', 'r') as f:
phi = mps.load_from_hdf5(config, f, './state/')
os.remove("tmp.h5")
assert (psi - phi).norm() < tol * psi.norm()