Initialization and creation operations#
Building MPS/MPO manually#
import numpy as np
import yastn
import yastn.tn.mps as mps
config_kwargs = {"backend": "np"}
The tensor making up MPS/MPO can be assigned manually, setting them one by one.
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.
def build_aklt_state_manually(config_kwargs, N=5, lvec=(1, 0), rvec=(0, 1)):
"""
Example for Spin-1 AKLT state of ``N`` sites.
Initialize MPS tensors by hand.
Allow inputing boundary vectors ``lvec`` an ``rvec`` for an MPS with OBC.
"""
# 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
#
config = yastn.make_config(sym='none', **config_kwargs)
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\), for different realizations of explicit symmetry: dense, Z2, or U1.
def mpo_nn_hopping_manually(config_kwargs, sym='U1', N=9, t=1, mu=1):
"""
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'.
"""
#
# 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.
config = yastn.make_config(sym=sym, fermionic=True, **config_kwargs)
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 mpo_hopping_Hterm(config_kwargs, sym='U1', J=[[1, 1, 1], [0, 1, 1], [0, 0, 1]]):
"""
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.
"""
ops = yastn.operators.SpinlessFermions(sym=sym, **config_kwargs)
c, cp, occ = ops.c(), ops.cp(), ops.n()
#
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=occ))
#
# 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=(cp, c)))
Hterms.append(mps.Hterm(amplitude=np.conj(J[m][n]),
positions=(n, m),
operators=(cp, 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(config_kwargs, sym='Z2', N=10, D_total=16, n=1):
"""
Generate random MPS of N sites, with bond dimension D_total,
tensors with symmetry sym and total charge n.
"""
#
# load local operators
#
ops = yastn.operators.SpinlessFermions(sym=sym, **config_kwargs)
#
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(config_kwargs, sym='U1', N=9, t=1, mu=1):
"""
Nearest-neighbor hopping Hamiltonian on N sites
with hopping amplitude t and chemical potential mu.
"""
ops = yastn.operators.SpinlessFermions(sym=sym, **config_kwargs)
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": [(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#
@pytest.mark.parametrize('sym', ['dense', 'Z3', 'U1'])
def test_save_load_mps_dict(config_kwargs, sym, tol=1e-12):
"""
Initialize random MPS and checks saving/loading to/from npy file.
"""
#
# generate random mps with 3-dimensional local spaces
#
ops = yastn.operators.Spin1(sym=sym, **config_kwargs)
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 pytest
@pytest.mark.parametrize('sym', ['dense', 'Z3', 'U1'])
def test_save_load_mps_hdf5(config_kwargs, sym, tol=1e-12):
"""
Initialize random MPS and checks saving/loading to/from HDF5 file.
"""
h5py = pytest.importorskip('h5py')
#
# generate random mps with 3-dimensional local spaces
#
ops = yastn.operators.Spin1(sym=sym, **config_kwargs)
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()