r"""This module contains a base class for a Matrix Product State (MPS).
An MPS looks roughly like this::
| -- B[0] -- B[1] -- B[2] -- ...
| | | |
We use the following label convention for the `B` (where arrows indicate `qconj`)::
| vL ->- B ->- vR
| |
| ^
| p
We store one 3-leg tensor `_B[i]` with labels ``'vL', 'vR', 'p'`` for each of the `L` sites
``0 <= i < L``.
Additionally, we store ``L+1`` singular value arrays `_S[ib]` on each bond ``0 <= ib <= L``,
independent of the boundary conditions.
``_S[ib]`` gives the singlur values on the bond ``i-1, i``.
However, be aware that e.g. :attr:`~tenpy.networks.mps.MPS.chi` returns only the dimensions of the
:attr:`~tenpy.networks.mps.MPS.nontrivial_bonds` depending on the boundary conditions.
The matrices and singular values always represent a normalized state
(i.e. ``np.linalg.norm(psi._S[ib]) == 1`` up to roundoff errors),
but (for finite MPS) we keep track of the norm in :attr:`~tenpy.networks.mps.MPS.norm`
(which is respected by :meth:`~tenpy.networks.mps.MPS.overlap`, ...).
Valid MPS boundary conditions (not to confuse with `bc_coupling` of
:class:`tenpy.models.model.CouplingModel`) are the following:
========== ===================================================================================
`bc` description
========== ===================================================================================
'finite' Finite MPS, ``G0 s1 G1 ... s{L-1} G{l-1}``. This is acchieved
by using a trivial left and right bond ``s[0] = s[-1] = np.array([1.])``.
'segment' Generalization of 'finite', describes an MPS embedded in left and right
environments. The left environment is described by ``chi[0]`` *orthonormal* states
which are weighted by the singular values ``s[0]``. Similar, ``s[L]`` weight some
right orthonormal states. You can think of the left and right states to be
generated by additional MPS, such that the overall structure is something like
``... s L s L [s0 G0 s1 G1 ... s{L-1} G{L-1} s{L}] R s R s R ...``
(where we save the part in the brackets ``[ ... ]`` ).
'infinite' infinite MPS (iMPS): we save a 'MPS unit cell' ``[s0 G0 s1 G1 ... s{L-1} G{L-1}]``
which is repeated periodically, identifying all indices modulo ``self.L``.
In particular, the last bond ``L`` is identified with ``0``.
(The MPS unit cell can differ from a lattice unit cell).
bond is identified with the first one.
========== ===================================================================================
An MPS can be in different 'canonical forms' (see [Vidal2004]_, [Schollwoeck2011]_).
To take care of the different canonical forms, algorithms should use functions like
:meth:`~tenpy.networks.mps.MPS.get_theta`, :meth:`~tenpy.networks.mps.MPS.get_B`
and :meth:`~tenpy.networks.mps.MPS.set_B` instead of accessing them directly,
as they return the `B` in the desired form (which can be chosen as an argument).
The values of the tuples for the form correspond to the exponent of the singular values
on the left and right.
To keep track of a "mixed" canonical form ``A A A s B B``, we save the tuples for each
site of the MPS in :attr:`MPS.form`.
======== ========== ==========================================================================
`form` tuple description
======== ========== ==========================================================================
``'B'`` (0, 1) right canonical: ``_B[i] = -- Gamma[i] -- s[i+1]--``
The default form, which algorithms asssume.
``'C'`` (0.5, 0.5) symmetric form: ``_B[i] = -- s[i]**0.5 -- Gamma[i] -- s[i+1]**0.5--``
``'A'`` (1, 0) left canonical: ``_B[i] = -- s[i] -- Gamma[i] --``.
``'G'`` (0, 0) Save only ``_B[i] = -- Gamma[i] --``.
``'Th'`` (1, 1) Form of a local wave function `theta` with singular value on both sides.
``psi.get_B(i, 'Th') is equivalent to ``psi.get_theta(i, n=1)``.
``None`` ``None`` General non-canoncial form.
Valid form for initialization, but you need to call
:meth:`~tenpy.networks.mps.MPS.canonical_form` (or similar)
before using algorithms.
======== ========== ==========================================================================
"""
# Copyright 2018-2019 TeNPy Developers, GNU GPLv3
import numpy as np
import warnings
import random
from functools import reduce
import scipy.sparse.linalg.eigen.arpack
from ..linalg import np_conserved as npc
from ..linalg import sparse
from .site import GroupedSite, group_sites
from ..tools.misc import to_iterable, argsort
from ..tools.math import lcm, speigs, entropy
from ..algorithms.truncation import TruncationError, svd_theta
__all__ = ['MPS', 'MPSEnvironment', 'TransferMatrix', 'build_initial_state']
[docs]class MPS:
r"""A Matrix Product State, finite (MPS) or infinite (iMPS).
Parameters
----------
sites : list of :class:`~tenpy.networks.site.Site`
Defines the local Hilbert space for each site.
Bs : list of :class:`~tenpy.linalg.np_conserved.Array`
The 'matrices' of the MPS. Labels are ``vL, vR, p`` (in any order).
SVs : list of 1D array
The singular values on *each* bond. Should always have length `L+1`.
Entries out of :attr:`nontrivial_bonds` are ignored.
bc : ``'finite' | 'segment' | 'infinite'``
Boundary conditions as described in the tabel of the module doc-string.
form : (list of) {``'B' | 'A' | 'C' | 'G' | 'Th' | None`` | tuple(float, float)}
The form of the stored 'matrices', see table in module doc-string.
A single choice holds for all of the entries.
Attributes
----------
L
chi
finite
nontrivial_bonds
sites : list of :class:`~tenpy.networks.site.Site`
Defines the local Hilbert space for each site.
bc : {'finite', 'segment', 'infinite'}
Boundary conditions as described in above table.
form : list of {``None`` | tuple(float, float)}
Describes the canonical form on each site.
``None`` means non-canonical form.
For ``form = (nuL, nuR)``, the stored ``_B[i]`` are
``s**form[0] -- Gamma -- s**form[1]`` (in Vidal's notation).
chinfo : :class:`~tenpy.linalg.np_conserved.ChargeInfo`
The nature of the charge.
dtype : type
The data type of the ``_B``.
norm : float
The norm of the state, i.e. ``sqrt(<psi|psi>)``.
Ignored for (normalized) :meth:`expectation_value`, but important for :meth:`overlap`.
grouped : int
Number of sites grouped together, see :meth:`group_sites`.
_B : list of :class:`npc.Array`
The 'matrices' of the MPS. Labels are ``vL, vR, p`` (in any order).
We recommend using :meth:`get_B` and :meth:`set_B`, which will take care of the different
canonical forms.
_S : list of (``None`` | 1D array)
The singular values on each virtual bond, length ``L+1``.
May be ``None`` if the MPS is not in canonical form.
Otherwise, ``_S[i]`` is to the left of ``_B[i]``.
We recommend using :meth:`get_SL`, :meth:`get_SR`, :meth:`set_SL`, :meth:`set_SR`, which
takes proper care of the boundary conditions.
_valid_forms : dict
Mapping for canonical forms to a tuple ``(nuL, nuR)`` indicating that
``self._Bs[i] = s[i]**nuL -- Gamma[i] -- s[i]**nuR`` is saved.
_valid_bc : tuple of str
Valid boundary conditions.
_transfermatrix_keep : int
How many states to keep at least when diagonalizing a :class:`TransferMatrix`.
Important if the state develops a near-degeneracy.
"""
# Canonical form conventions: the saved B = s**nu[0]--Gamma--s**nu[1].
# For the canonical forms, ``nu[0] + nu[1] = 1``
_valid_forms = {
'A': (1., 0.),
'C': (0.5, 0.5),
'B': (0., 1.),
'G': (0., 0.), # like Vidal's `Gamma`.
'Th': (1., 1.),
None: None, # means 'not in any canonical form'
}
# valid boundary conditions. Don't overwrite this!
_valid_bc = ('finite', 'segment', 'infinite')
# the "physical" labels for each B
_p_label = ['p']
_p_label_star = ['p*']
# All labels of each tensor in _B (order is used!)
_B_labels = ['vL', 'p', 'vR']
def __init__(self, sites, Bs, SVs, bc='finite', form='B', norm=1.):
self.sites = list(sites)
self.chinfo = self.sites[0].leg.chinfo
self.dtype = dtype = np.find_common_type([B.dtype for B in Bs], [])
self.form = self._parse_form(form)
self.bc = bc # one of ``'finite', 'periodic', 'segment'``.
self.norm = norm
self.grouped = 1
# make copies of Bs and SVs
self._B = [B.astype(dtype, copy=True).itranspose(self._B_labels) for B in Bs]
self._S = [None] * (self.L + 1)
for i in range(self.L + 1)[self.nontrivial_bonds]:
self._S[i] = np.array(SVs[i], dtype=np.float)
if self.bc == 'infinite':
self._S[-1] = self._S[0]
elif self.bc == 'finite':
self._S[0] = self._S[-1] = np.ones([1])
self._transfermatrix_keep = 1
self.test_sanity()
[docs] def test_sanity(self):
"""Sanity check, raises ValueErrors, if something is wrong."""
if self.bc not in self._valid_bc:
raise ValueError("invalid boundary condition: " + repr(self.bc))
if len(self._B) != self.L:
raise ValueError("wrong len of self._B")
if len(self._S) != self.L + 1:
raise ValueError("wrong len of self._S")
for i, B in enumerate(self._B):
if B.get_leg_labels() != self._B_labels:
raise ValueError("B has wrong labels {0!r}, expected {1!r}".format(
B.get_leg_labels(), self._B_labels))
if self._S[i].shape[-1] != B.get_leg('vL').ind_len or \
self._S[i+1].shape[0] != B.get_leg('vR').ind_len:
raise ValueError("shape of B incompatible with len of singular values")
if not self.finite or i + 1 < self.L:
B2 = self._B[(i + 1) % self.L]
B.get_leg('vR').test_contractible(B2.get_leg('vL'))
if self.bc == 'finite':
if len(self._S[0]) != 1 or len(self._S[-1]) != 1:
raise ValueError("non-trivial outer bonds for finite MPS")
elif self.bc == 'infinite':
if np.any(self._S[self.L] != self._S[0]):
raise ValueError("iMPS with S[0] != S[L]")
assert len(self.form) == self.L
for f in self.form:
if f is not None:
assert isinstance(f, tuple)
assert len(f) == 2
[docs] @classmethod
def from_product_state(cls,
sites,
p_state,
bc='finite',
dtype=np.float64,
permute=True,
form='B',
chargeL=None):
"""Construct a matrix product state from a given product state.
Parameters
----------
sites : list of :class:`~tenpy.networks.site.Site`
The sites defining the local Hilbert space.
p_state : iterable of {int | str | 1D array}
Defines the product state to be represented.
If ``p_state[i]`` is str, then site ``i`` is in state
``self.sites[i].state_labels(p_state[i])``.
If ``p_state[i]`` is int, then site ``i`` is in state ``p_state[i]``.
If ``p_state[i]`` is an array, then site ``i`` wavefunction is ``p_state[i]``.
bc : {'infinite', 'finite', 'segmemt'}
MPS boundary conditions. See docstring of :class:`MPS`.
dtype : type or string
The data type of the array entries.
permute : bool
The :class:`~tenpy.networks.Site` might permute the local basis states if charge
conservation gets enabled.
If `permute` is True (default), we permute the given `p_state` locally according to
each site's :attr:`~tenpy.networks.Site.perm`.
The `p_state` argument should then always be given as if `conserve=None` in the Site.
form : (list of) {``'B' | 'A' | 'C' | 'G' | None`` | tuple(float, float)}
Defines the canonical form. See module doc-string.
A single choice holds for all of the entries.
chargeL : charges
Leg charge at bond 0, which are purely conventional.
Returns
-------
product_mps : :class:`MPS`
An MPS representing the specified product state.
Examples
--------
Example to get a Neel state for a :class:`~tenpy.models.tf_ising.TIChain`:
>>> M = TFIChain({'L': 10})
>>> p_state = ["up", "down"] * (L//2) # repeats entries L/2 times
>>> psi = MPS.from_product_state(M.lat.mps_sites(), p_state, bc=M.lat.bc_MPS)
The meaning of the labels ``"up","down"`` is defined by the :class:`~tenpy.networks.Site`,
in this example a :class:`~tenpy.networks.site.SpinHalfSite`.
Extending the example, we can replace the spin in the center with one with arbitrary
angles ``theta, phi`` in the bloch sphere:
>>> M = TFIChain({'L': 8, 'conserve': None})
>>> p_state = ["up", "down"] * (L//2) # repeats entries L/2 times
>>> bloch_sphere_state = np.array([np.cos(theta/2), np.exp(1.j*phi)*np.sin(theta/2)])
>>> p_state[L//2] = bloch_sphere_state # replace one spin in center
>>> psi = MPS.from_product_state(M.lat.mps_sites(), p_state, bc=M.lat.bc_MPS, dtype=np.complex)
Note that for the more general :class:`~tenpy.models.spins.SpinChain`,
the order of the two entries for the ``bloch_sphere_state`` would be *exactly the opposite*
(when we keep the the north-pole of the bloch sphere being the up-state).
The reason is that the `SpinChain` uses the general :class:`~tenpy.networks.site.SpinSite`,
where the states are orderd ascending from ``'down'`` to ``'up'``.
The :class:`~tenpy.networks.site.SpinHalfSite` on the other hand uses the order
``'up', 'down'`` where that the Pauli matrices look as usual.
Moreover, note that you can not write this bloch state (for ``theta != 0, pi``) when
conserving symmetries, as the two physical basis states correspond to different symmetry
sectors.
"""
sites = list(sites)
L = len(sites)
p_state = list(p_state)
if len(p_state) != L:
raise ValueError("Length of p_state does not match number of sites.")
ci = sites[0].leg.chinfo
Bs = []
chargeL = ci.make_valid(chargeL) # sets to zero if `None`
legL = npc.LegCharge.from_qflat(ci, [chargeL]) # (no need to bunch)
for p_st, site in zip(p_state, sites):
perm = permute
if isinstance(p_st, str):
p_st = site.state_labels[p_st] # translate labels into "int"
perm = False
try:
iter(p_st)
except TypeError:
# just an int for p_st
B = np.zeros((site.dim, 1, 1), dtype)
B[p_st, 0, 0] = 1.0
else: # iter works
if len(p_st) != site.dim:
raise ValueError("p_state incompatible with local dim:" + repr(p_st))
B = np.array(p_st, dtype).reshape((site.dim, 1, 1))
if perm:
B = B[site.perm, :, :]
Bs.append(B)
SVs = [[1.]] * (L + 1)
return cls.from_Bflat(sites, Bs, SVs, bc, dtype, False, form, legL)
[docs] @classmethod
def from_Bflat(cls,
sites,
Bflat,
SVs=None,
bc='finite',
dtype=None,
permute=True,
form='B',
legL=None):
"""Construct a matrix product state from a set of numpy arrays `Bflat` and singular vals.
Parameters
----------
sites : list of :class:`~tenpy.networks.site.Site`
The sites defining the local Hilbert space.
Bflat : iterable of numpy ndarrays
The matrix defining the MPS on each site, with legs ``'p', 'vL', 'vR'``
(physical, virtual left/right).
SVs : list of 1D array | ``None``
The singular values on *each* bond. Should always have length `L+1`.
By default (``None``), set all singular values to the same value.
Entries out of :attr:`nontrivial_bonds` are ignored.
bc : {'infinite', 'finite', 'segmemt'}
MPS boundary conditions. See docstring of :class:`MPS`.
dtype : type or string
The data type of the array entries. Defaults to the common dtype of `Bflat`.
permute : bool
The :class:`~tenpy.networks.Site` might permute the local basis states if charge
conservation gets enabled.
If `permute` is True (default), we permute the given `Bflat` locally according to
each site's :attr:`~tenpy.networks.Site.perm`.
The `p_state` argument should then always be given as if `conserve=None` in the Site.
form : (list of) {``'B' | 'A' | 'C' | 'G' | None`` | tuple(float, float)}
Defines the canonical form of `Bflat`. See module doc-string.
A single choice holds for all of the entries.
leg_L : LegCharge | ``None``
Leg charges at bond 0, which are purely conventional.
If ``None``, use trivial charges.
Returns
-------
mps : :class:`MPS`
An MPS with the matrices `Bflat` converted to npc arrays.
"""
sites = list(sites)
L = len(sites)
Bflat = list(Bflat)
if len(Bflat) != L:
raise ValueError("Length of Bflat does not match number of sites.")
ci = sites[0].leg.chinfo
if legL is None:
legL = npc.LegCharge.from_qflat(ci, [ci.make_valid(None)] * Bflat[0].shape[1])
legL = legL.bunch()[1]
if SVs is None:
SVs = [np.ones(B.shape[1]) / np.sqrt(B.shape[1]) for B in Bflat]
SVs.append(np.ones(Bflat[-1].shape[2]) / np.sqrt(Bflat[-1].shape[2]))
Bs = []
if dtype is None:
dtype = np.dtype(np.common_type(*Bflat))
for i, site in enumerate(sites):
B = np.array(Bflat[i], dtype)
if permute:
B = B[site.perm, :, :]
# calculate the LegCharge of the right leg
legs = [site.leg, legL, None] # other legs are known
legs = npc.detect_legcharge(B, ci, legs, None, qconj=-1)
B = npc.Array.from_ndarray(B, legs, dtype)
B.iset_leg_labels(['p', 'vL', 'vR'])
Bs.append(B)
legL = legs[-1].conj() # prepare for next `i`
if bc == 'infinite':
# for an iMPS, the last leg has to match the first one.
# so we need to gauge `qtotal` of the last `B` such that the right leg matches.
chdiff = Bs[-1].get_leg('vR').charges[0] - Bs[0].get_leg('vL').charges[0]
Bs[-1] = Bs[-1].gauge_total_charge('vR', ci.make_valid(chdiff))
return cls(sites, Bs, SVs, form=form, bc=bc)
[docs] @classmethod
def from_full(cls,
sites,
psi,
form=None,
cutoff=1.e-16,
normalize=True,
bc='finite',
outer_S=None):
"""Construct an MPS from a single tensor `psi` with one leg per physical site.
Performs a sequence of SVDs of psi to split off the `B` matrices and obtain the singular
values, the result will be in canonical form.
Obviously, this is only well-defined for `finite` or `segment` boundary conditions.
Parameters
----------
sites : list of :class:`~tenpy.networks.site.Site`
The sites defining the local Hilbert space.
psi : :class:`~tenpy.linalg.np_conserved.Array`
The full wave function to be represented as an MPS.
Should have labels ``'p0', 'p1', ..., 'p{L-1}'``.
Additionally, it may have (or must have for 'segment' `bc`) the legs ``'vL', 'vR'``,
which are trivial for 'finite' `bc`.
form : ``'B' | 'A' | 'C' | 'G' | None``
The canonical form of the resulting MPS, see module doc-string.
``None`` defaults to 'A' form on the first site and 'B' form on all following sites.
cutoff : float
Cutoff of singular values used in the SVDs.
normalize : bool
Whether the resulting MPS should have 'norm' 1.
bc : 'finite' | 'segment'
Boundary conditions.
outer_S : None | (array, array)
For 'semgent' `bc` the singular values on the left and right of the considered segment,
`None` for 'finite' boundary conditions.
Returns
-------
psi_mps : :class:`MPS`
MPS representation of `psi`, in canonical form and possibly normalized.
"""
if form is not None and form not in ['B', 'A', 'C', 'G']:
raise ValueError("Invalid form: " + repr(form))
if bc != 'finite' and bc != 'segment':
raise ValueError("Wrong boundary conditions: " + repr(bc))
# perform SVDs to bring it into 'B' form, afterwards change the form.
L = len(sites)
assert (L >= 2)
B_list = [None] * L
S_list = [None] * (L + 1)
norm = 1. if normalize else npc.norm(psi)
if not psi.has_label('vL'):
psi = psi.add_trivial_leg(0, label='vL', qconj=+1)
elif bc == 'finite' and psi.get_leg('vL').ind_len != 1:
raise ValueError("non-trivial left leg for 'finite' bc!")
if not psi.has_label('vR'):
psi = psi.add_trivial_leg(len(psi.get_leg_labels()), label='vR', qconj=-1)
elif bc == 'finite' and psi.get_leg('vR').ind_len != 1:
raise ValueError("non-trivial left leg for 'finite' bc!")
labels = ['vL'] + ['p' + str(i) for i in range(L)] + ['vR']
psi.itranspose(labels)
# combine legs from left
for i in range(0, L - 1):
psi = psi.combine_legs([0, 1]) # combines the legs until `i`
# now psi has only three legs: ``'(((vL.p0).p1)...p{L-2})', 'p{L-1}', 'vR'``
for i in range(L - 1, 0, -1):
# split off B[i]
psi = psi.combine_legs([labels[i + 1], 'vR'])
psi, S, B = npc.svd(psi, inner_labels=['vR', 'vL'], cutoff=cutoff)
S /= np.linalg.norm(S) # normalize
if i > 1:
psi.iscale_axis(S, 1)
B_list[i] = B.split_legs(1).replace_label(labels[i + 1], 'p')
S_list[i] = S
psi = psi.split_legs(0)
# psi is now the first `B` in 'A' form
B_list[0] = psi.replace_label(labels[1], 'p')
B_form = ['A'] + ['B'] * (L - 1)
if bc == 'finite':
S_list[0] = S_list[-1] = np.ones([1], dtype=np.float)
elif outer_S is not None:
S_list[0], S_list[-1] = outer_S
res = cls(sites, B_list, S_list, bc=bc, form=B_form, norm=norm)
if form is not None:
res.convert_form(form)
return res
[docs] @classmethod
def from_singlets(cls,
site,
L,
pairs,
up='up',
down='down',
lonely=[],
lonely_state='up',
bc='finite'):
"""Create an MPS of entangled singlets.
Parameters
----------
site : :class:`~tenpy.networks.site.Site`
The `site` defining the local Hilbert space, taken uniformly for all sites.
L : int
The number of sites.
pairs : list of (int, int)
Pairs of sites to be entangled; the returned MPS will have a singlet
for each pair in `pairs`.
up, down : int | str
A singlet is defined as ``(|up down> - |down up>)/2**0.5``,
``up`` and ``down`` give state indices or labels defined on the corresponding site.
lonely : list of int
Sites which are not included into a singlet pair.
lonely_state : int | str
The state for the lonely sites.
bc : {'infinite', 'finite', 'segmemt'}
MPS boundary conditions. See docstring of :class:`MPS`.
Returns
-------
singlet_mps : :class:`MPS`
An MPS representing singlets on the specified pairs of sites.
"""
# sort each pair s.t. i < j
pairs = [((i, j) if i < j else (j, i)) for (i, j) in pairs]
# sort by smaller site of the pair
pairs.sort(key=lambda x: x[0])
pairs.append((L, L))
lonely = sorted(lonely) + [L]
# generate building block tensors
up = site.state_index(up)
down = site.state_index(down)
lonely_state = site.state_index(lonely_state)
mask = np.zeros(site.dim, dtype=np.bool_)
mask[up] = mask[down] = True
Open = npc.diag(1., site.leg)[:, mask]
Close = np.zeros([site.dim, site.dim], dtype=np.float_)
Close[up, down] = 1.
Close[down, up] = -1.
Close = npc.Array.from_ndarray(Close, [site.leg, site.leg]) # no conj() !
Close = Close[mask, :]
Id = npc.eye_like(Close, 0)
Lonely = np.zeros(site.dim, dtype=np.float_)
Lonely[lonely_state] = 1
Lonely = npc.Array.from_ndarray(Lonely, [site.leg])
Bs = []
Ss = [np.ones(1)]
forms = []
open_singlets = [] # the k-th open singlet should be closed at site open_singlets[k]
Ts = [] # the tensors on the current site
labels_L = []
for i in range(L):
labels_R = labels_L[:]
next_Ts = Ts[:]
if i == pairs[0][0]: # open a new singlet
j = pairs[0][1]
lbl = 's{0:d}-{1:d}'.format(i, j)
pairs.pop(0)
open_singlets.append(j)
next_Ts.append(Id.copy().iset_leg_labels([lbl + 'L', lbl]))
Open.iset_leg_labels(['p', lbl])
Ts.append(Open.copy(deep=False))
labels_R.append(lbl)
forms.append('A')
elif i == lonely[0]: # just a lonely state
Ts.append(Lonely)
lonely.pop(0)
forms.append('B')
else: # close a singlet
k = open_singlets.index(i)
Close.iset_leg_labels([labels_L[k] + 'L', 'p'])
Ts[k] = Close
next_Ts.pop(k)
open_singlets.pop(k)
labels_R.pop(k)
forms.append('B')
# generate `B` from `Ts`
B = reduce(npc.outer, Ts)
labels_L = [lbl_ + 'L' for lbl_ in labels_L]
if len(labels_L) > 0 and len(labels_R) > 0:
B = B.combine_legs([labels_L, labels_R], new_axes=[0, 2], qconj=[+1, -1])
B.iset_leg_labels(['vL', 'p', 'vR'])
elif len(labels_L) == 0 and len(labels_R) == 0:
B = B.add_trivial_leg(0, label='vL', qconj=+1)
B = B.add_trivial_leg(2, label='vR', qconj=-1)
B.iset_leg_labels(['vL', 'p', 'vR'])
elif len(labels_L) == 0:
B = B.combine_legs([labels_R], new_axes=[1], qconj=[-1])
B.iset_leg_labels(['p', 'vR'])
B = B.add_trivial_leg(0, label='vL', qconj=+1)
else: # len(labels_R) == 0
B = B.combine_legs([labels_L], new_axes=[0], qconj=[+1])
B.iset_leg_labels(['vL', 'p'])
B = B.add_trivial_leg(2, label='vR', qconj=-1)
Bs.append(B)
N = 2**len(labels_R)
Ss.append(np.ones(N) / (N**0.5))
Ts = next_Ts
labels_L = labels_R
return cls([site] * L, Bs, Ss, bc=bc, form=forms)
[docs] def copy(self):
"""Returns a copy of `self`.
The copy still shares the sites, chinfo, and LegCharges of the B tensors, but the values of
B and S are deeply copied.
"""
# __init__ makes deep copies of B, S
cp = MPS(self.sites, self._B, self._S, self.bc, self.form, self.norm)
cp.grouped = self.grouped
cp._transfermatrix_keep = self._transfermatrix_keep
return cp
@property
def L(self):
"""Number of physical sites; for an iMPS the len of the MPS unit cell."""
return len(self.sites)
@property
def dim(self):
"""List of local physical dimensions."""
return [site.dim for site in self.sites]
@property
def finite(self):
"""Distinguish MPS vs iMPS.
True for an MPS (``bc='finite', 'segment'``), False for an iMPS (``bc='infinite'``).
"""
assert (self.bc in self._valid_bc)
return self.bc != 'infinite'
@property
def chi(self):
"""Dimensions of the (nontrivial) virtual bonds."""
# s.shape[0] == len(s) for 1D numpy array, but works also for a 2D npc Array.
return [min(s.shape) for s in self._S[self.nontrivial_bonds]]
@property
def nontrivial_bonds(self):
"""Slice of the non-trivial bond indices, depending on ``self.bc``."""
if self.bc == 'finite':
return slice(1, self.L)
elif self.bc == 'segment':
return slice(0, self.L + 1)
elif self.bc == 'infinite':
return slice(0, self.L)
[docs] def get_B(self, i, form='B', copy=False, cutoff=1.e-16, label_p=None):
"""Return (view of) `B` at site `i` in canonical form.
Parameters
----------
i : int
Index choosing the site.
form : ``'B' | 'A' | 'C' | 'G' | 'Th' | None`` | tuple(float, float)
The (canonical) form of the returned B.
For ``None``, return the matrix in whatever form it is.
If any of the tuple entry is None, also don't scale on the corresponding axis.
copy : bool
Whether to return a copy even if `form` matches the current form.
cutoff : float
During DMRG with a mixer, `S` may be a matrix for which we need the inverse.
This is calculated as the Penrose pseudo-inverse, which uses a cutoff for the
singular values.
label_p : None | str
Ignored by default (``None``).
Otherwise replace the physical label ``'p'`` with ``'p'+label_p'``.
(For derived classes with more than one "physical" leg, replace all the physical leg
labels accordingly.)
Returns
-------
B : :class:`~tenpy.linalg.np_conserved.Array`
The MPS 'matrix' `B` at site `i` with leg labels ``'vL', 'p', 'vR'``.
May be a view of the matrix (if ``copy=False``),
or a copy (if the form changed or ``copy=True``).
Raises
------
ValueError : if self is not in canoncial form and `form` is not None.
"""
i = self._to_valid_index(i)
new_form = self._to_valid_form(form)
old_form = self.form[i]
B = self._B[i]
if copy:
B = B.copy()
if new_form is not None and old_form != new_form:
if old_form is None:
raise ValueError("can't convert form of non-canonical state!")
if new_form[0] is not None and new_form[0] - old_form[0] != 0.:
B = self._scale_axis_B(B, self.get_SL(i), new_form[0] - old_form[0], 'vL', cutoff)
if new_form[1] is not None and new_form[1] - old_form[1] != 0.:
B = self._scale_axis_B(B, self.get_SR(i), new_form[1] - old_form[1], 'vR', cutoff)
if label_p is not None:
B = self._replace_p_label(B, label_p)
return B
[docs] def set_B(self, i, B, form='B'):
"""Set `B` at site `i`.
Parameters
----------
i : int
Index choosing the site.
B : :class:`~tenpy.linalg.np_conserved.Array`
The 'matrix' at site `i`. No copy is made!
Should have leg labels ``'vL', 'p', 'vR'`` (not necessarily in that order).
form : ``'B' | 'A' | 'C' | 'G' | 'Th' | None`` | tuple(float, float)
The (canonical) form of the `B` to set.
``None`` stands for non-canonical form.
"""
i = self._to_valid_index(i)
self.form[i] = self._to_valid_form(form)
self.dtype = np.find_common_type([self.dtype, B.dtype], [])
self._B[i] = B.itranspose(self._B_labels)
[docs] def get_SL(self, i):
"""Return singular values on the left of site `i`"""
i = self._to_valid_index(i)
return self._S[i]
[docs] def get_SR(self, i):
"""Return singular values on the right of site `i`"""
i = self._to_valid_index(i)
return self._S[i + 1]
[docs] def set_SL(self, i, S):
"""Set singular values on the left of site `i`"""
i = self._to_valid_index(i)
self._S[i] = S
if not self.finite and i == 0:
self._S[self.L] = S
[docs] def set_SR(self, i, S):
"""Set singular values on the right of site `i`"""
i = self._to_valid_index(i)
self._S[i + 1] = S
if not self.finite and i == self.L - 1:
self._S[0] = S
[docs] def get_op(self, op_list, i):
"""Given a list of operators, select the one corresponding to site `i`.
Parameters
----------
op_list : (list of) {str | npc.array}
List of operators from which we choose. We assume that ``op_list[j]`` acts on site
``j``. If the length is shorter than `L`, we repeat it periodically.
Strings are translated using :meth:`~tenpy.networks.site.Site.get_op` of site `i`.
i : int
Index of the site on which the operator acts.
Returns
-------
op : npc.array
One of the entries in `op_list`, not copied.
"""
i = self._to_valid_index(i)
op = op_list[i % len(op_list)]
if (isinstance(op, str)):
op = self.sites[i].get_op(op)
return op
[docs] def get_theta(self, i, n=2, cutoff=1.e-16, formL=1., formR=1.):
"""Calculates the `n`-site wavefunction on ``sites[i:i+n]``.
Parameters
----------
i : int
Site index.
n : int
Number of sites. The result lives on ``sites[i:i+n]``.
cutoff : float
During DMRG with a mixer, `S` may be a matrix for which we need the inverse.
This is calculated as the Penrose pseudo-inverse, which uses a cutoff for the
singular values.
formL : float
Exponent for the singular values to the left.
formR : float
Exponent for the singular values to the right.
Returns
-------
theta : :class:`~tenpy.linalg.np_conserved.Array`
The n-site wave function with leg labels ``vL, p0, p1, .... p{n-1}, vR``.
In Vidal's notation (with s=lambda, G=Gamma):
``theta = s**form_L G_i s G_{i+1} s ... G_{i+n-1} s**form_R``.
"""
i = self._to_valid_index(i)
for j in range(i, i + n):
if self.form[j % self.L] is None:
raise ValueError("can't calculate theta for non-canonical form")
if n == 1:
return self.get_B(i, (1., 1.), True, cutoff, '0')
# n >= 2: contract some B's
theta = self.get_B(i, (formL, None), False, cutoff, '0') # right form as stored
_, old_fR = self.form[i]
for k in range(1, n): # non-empty range
j = self._to_valid_index(i + k)
new_fR = None if k + 1 < n else formR # right form as stored, except for last B
B = self.get_B(j, (1. - old_fR, new_fR), False, cutoff, str(k))
_, old_fR = self.form[j]
theta = npc.tensordot(theta, B, axes=['vR', 'vL'])
return theta
[docs] def increase_L(self, new_L=None):
"""Modify `self` inplace to enlarge the unit cell.
For an infinite MPS, we have unit cells.
Parameters
----------
new_L : int
New number of sites. Defaults to twice the number of current sites.
"""
old_L = self.L
if new_L is None:
new_L = 2 * old_L
if new_L % old_L:
raise ValueError("new_L = {0:d} not a multiple of old L={1:d}".format(new_L, old_L))
factor = new_L // old_L
if factor <= 1:
raise ValueError("can't shrink!")
if self.bc == 'segment':
raise ValueError("can't enlarge segment MPS")
self.sites = factor * self.sites
self._B = factor * self._B
self._S = factor * self._S[:-1] + [self._S[-1]]
self.form = factor * self.form
[docs] def group_sites(self, n=2, grouped_sites=None):
"""Modify `self` inplace to group sites.
Group each `n` sites together using the :class:`~tenpy.networks.site.GroupedSite`.
This might allow to do TEBD with a Trotter decomposition,
or help the convergence of DMRG (in case of too long range interactions).
Parameters
----------
n : int
Number of sites to be grouped together.
grouped_sites : None | list of :class:`~tenpy.networks.site.GroupedSite`
The sites grouped together.
See also
--------
group_split : Reverts the grouping.
"""
self.convert_form('B')
if grouped_sites is None:
grouped_sites = group_sites(self.sites, n, charges='same')
else:
assert grouped_sites[0].n_sites == n
Bs = []
Ss = []
i = 0
B_form = self._valid_forms['B']
for gs in grouped_sites:
n_sites = gs.n_sites
new_B = self.get_theta(i, gs.n_sites, formL=B_form[0], formR=B_form[1])
comb_legs = [[lbl + str(k) for k in range(n_sites)] for lbl in self._p_label]
# comb_legs = [['p0', 'p1', ... ]] for usual MPS
axes = list(range(1, 1 + len(self._p_label))) # [1]
new_B = new_B.combine_legs(comb_legs, new_axes=axes, qconj=[+1] * len(axes))
new_B.legs[1].test_equal(gs.leg) # test legcharge compatibility
Bs.append(new_B.iset_leg_labels(self._B_labels)) # ['vL', 'p', 'vR']
Ss.append(self._S[i])
i += n_sites
Ss.append(self._S[-1]) # right-most singular values: need L+1 entries
self._B = Bs
self._S = Ss
self.sites = grouped_sites
self.form = [B_form] * len(grouped_sites)
self.grouped = self.grouped * n
[docs] def group_split(self, trunc_par=None):
"""Modify `self` inplace to split previously grouped sites.
Parameters
----------
trunc_par : dict
Parameters for truncation, see :func:`~tenpy.algorithms.truncation.truncate`.
Defaults to ``{'chi_max': max(self.chi)}``.
Returns
-------
trunc_err : :class:`~tenpy.algorithms.truncation.TruncationError`
The error introduced by the truncation for the splitting.
See also
--------
group_sites : Should have been used before to combine sites.
"""
if trunc_par is None:
trunc_par = {}
self.convert_form('B')
if self.L > 1:
trunc_par.setdefault('chi_max', max(self.chi))
n0 = self.sites[0].n_sites
sites = []
Bs = []
Ss = []
trunc_err = TruncationError()
for i, gs in enumerate(self.sites):
sites.extend(gs.sites)
n = gs.n_sites
Ss_new = []
Bs_new = []
B_gr = self.get_B(i, form='B').transpose(self._B_labels) # vL, p, vR
B_gr.idrop_labels(self._p_label) # avoid warning: split label not called '(...)'
n_p_label = len(self._p_label)
split_legs = list(range(1, 1 + n_p_label))
transp = [i for k in range(n) for i in range(1 + k, 1 + n * n_p_label, n)]
transp = ['vL'] + transp + ['vR']
B_gr = B_gr.split_legs(split_legs).itranspose(transp)
theta = self.get_theta(i, n=1)
theta.idrop_labels(self._get_p_label('0')) # avoid warning
theta = theta.split_legs(split_legs).itranspose(transp)
# for usual MPS, B_gr and theta have legs vL p0 p1 ... p{n-1} vR
# for MPS with legs p, q, they have legs vL p0 q0 p1 q1 ... q{-n-1} vR
combine = [list(range(B_gr.rank - n_p_label - 1)), list(range(-n_p_label - 1, 0, +1))]
# combine = [[0, 1, .... {n-1}], [-2, -1]] for usual MPS
axes_contr = [combine[1], list(range(1, 2 + n_p_label))]
for j in range(n - 1, 0, -1):
# split off the right-most physical leg and vR from theta
# theta: vL p0 ... pj vR
theta = theta.combine_legs(combine, qconj=[+1, -1])
U, S, V, err, _ = svd_theta(theta, trunc_par, inner_labels=['vR', 'vL'])
Ss_new.append(S)
trunc_err += err
theta = U.split_legs(0) # vL p0 ... pj-1 vR
for _ in range(n_p_label):
combine[0].pop()
B = V.split_legs(1).iset_leg_labels(self._B_labels) # vL p vR
B_gr = npc.tensordot(B_gr, B.conj(), axes=axes_contr) # vL p0 ... pj-1 vR
Bs_new.append(B)
Bs_new.append(B_gr.iset_leg_labels(self._B_labels)) # inversion free :)
Ss_new.append(self.get_SL(i))
Bs.extend(Bs_new[::-1])
Ss.extend(Ss_new[::-1])
Ss.append(self._S[-1])
self.sites = sites
self._B = Bs
self._S = Ss
self.grouped = self.grouped // n0
self.form = [self._valid_forms['B']] * len(sites)
self.test_sanity()
return trunc_err
[docs] def get_grouped_mps(self, blocklen):
r"""contract blocklen subsequent tensors into a single one and return result as a new MPS.
blocklen = number of subsequent sites to be combined.
Returns
-------
new MPS object with bunched sites.
"""
groupedMPS = self.copy()
groupedMPS.group_sites(n=blocklen)
return groupedMPS
[docs] def get_total_charge(self, only_physical_legs=False):
"""Calculate and return the `qtotal` of the whole MPS (when contracted).
Parameters
----------
only_physical_legs : bool
For ``'finite'`` boundary conditions, the total charge can be gauged away
by changing the LegCharge of the trivial legs on the left and right of the MPS.
This option allows to project out the trivial legs to get the actual "physical"
total charge.
Returns
-------
qtotal : charges
The sum of the `qtotal` of the individual `B` tensors.
"""
qtotal = np.sum([B.qtotal for B in self._B], axis=0)
if only_physical_legs:
if self.bc != 'finite':
raise ValueError("`only_physical_legs` not supported for bc=" + repr(self.bc))
qtotal += self._B[0].get_leg('vL').get_charge(0)
qtotal += self._B[-1].get_leg('vR').get_charge(0) # takes qconj into account
return self.chinfo.make_valid(qtotal)
[docs] def gauge_total_charge(self, qtotal=None, vL_leg=None, vR_leg=None):
"""Gauge the legcharges of the virtual bonds such that the MPS has a total `qtotal`.
Parameters
----------
qtotal : (list of) charges
If a single set of charges is given, it is the desired total charge of the MPS
(which :meth:`get_total_charge` will return afterwards).
By default (``None``), use 0 charges, unless vL_leg and vR_leg are specified, in which
case we adjust the total charge to match these legs.
vL_leg : None | LegCharge
Desired new virtual leg on the very left. Needs to have the same block strucuture as
current leg, but can have shifted charge entries.
vR_leg : None | LegCharge
Desired new virtual leg on the very right. Needs to have the same block strucuture as
current leg, but can have shifted charge entries.
Should be `vL_leg.conj()` for infinite MPS, if `qtotal` is not given.
"""
if self.chinfo.qnumber == 0:
return
if vL_leg is not None:
vL_chdiff = vL_leg.get_charge(0) - self._B[0].get_leg('vL').get_charge(0)
if vR_leg is not None:
vR_chdiff = vR_leg.get_charge(0) - self._B[-1].get_leg('vR').get_charge(0)
if qtotal is None:
if vL_leg is not None and vR_leg is not None:
qtotal = self.get_total_charge() + vL_chdiff + vR_chdiff
qtotal = self.chinfo.make_valid(qtotal)
if qtotal.ndim == 1:
qtotal_factor = np.array([0] * (self.L - 1) + [1], npc.QTYPE)
qtotal = qtotal_factor[:, np.newaxis] * qtotal[np.newaxis, :]
if qtotal.shape != (self.L, self.chinfo.qnumber):
raise ValueError("wrong shape of `qtotal`")
if vL_leg is not None:
B = self._B[0]
if np.any(vL_chdiff != 0):
# adjust left leg
self._B[0] = B.gauge_total_charge('vL', B.qtotal + vL_chdiff, vL_leg.qconj)
B.get_leg('vL').test_equal(vL_leg)
for i in range(self.L):
B = self._B[i]
desired_qtotal = qtotal[i]
chdiff = B.qtotal - desired_qtotal
if np.any(chdiff != 0):
self._B[i] = B.gauge_total_charge('vR', desired_qtotal)
if i + 1 != self.L: # this 'vR' is contracted with the 'vL' of the next B
# so we need to adjust the next B as well
nextB = self._B[i + 1]
self._B[i + 1] = nextB.gauge_total_charge('vL', nextB.qtotal + chdiff)
self._B[i].get_leg('vR').test_contractible(self._B[i + 1].get_leg('vL'))
# just to check
assert np.all(self.get_total_charge() == self.chinfo.make_valid(np.sum(qtotal, 0)))
if vR_leg is not None:
# check that the charges match
self._B[-1].get_leg('vR').test_equal(vR_leg)
if self.bc == 'infinite':
self._B[0].get_leg('vL').test_contractable(self._B[-1].get_leg('vR'))
# done
[docs] def entanglement_entropy(self, n=1, bonds=None, for_matrix_S=False):
r"""Calculate the (half-chain) entanglement entropy for all nontrivial bonds.
Consider a bipartition of the sytem into :math:`A = \{ j: j <= i_b \}` and
:math:`B = \{ j: j > i_b\}` and the reduced density matrix :math:`\rho_A = tr_B(\rho)`.
The von-Neumann entanglement entropy is defined as
:math:`S(A, n=1) = -tr(\rho_A \log(\rho_A)) = S(B, n=1)`.
The generalization for ``n != 1, n>0`` are the Renyi entropies:
:math:`S(A, n) = \frac{1}{1-n} \log(tr(\rho_A^2)) = S(B, n=1)`
This function calculates the entropy for a cut at different bonds `i`, for which the
the eigenvalues of the reduced density matrix :math:`\rho_A` and :math:`\rho_B` is given
by the squared schmidt values `S` of the bond.
Parameters
----------
n : int/float
Selects which entropy to calculate;
`n=1` (default) is the ususal von-Neumann entanglement entropy.
bonds : ``None`` | (iterable of) int
Selects the bonds at which the entropy should be calculated.
``None`` defaults to ``range(0, L+1)[self.nontrivial_bonds]``.
for_matrix_S : bool
Switch calculate the entanglement entropy even if the `_S` are matrices.
Since :math:`O(\chi^3)` is expensive compared to the ususal :math:`O(\chi)`,
we raise an error by default.
Returns
-------
entropies : 1D ndarray
Entanglement entropies for half-cuts.
`entropies[j]` contains the entropy for a cut at bond ``bonds[j]``
(i.e. left to site ``bonds[j]``).
"""
if bonds is None:
nt = self.nontrivial_bonds
bonds = range(nt.start, nt.stop)
res = []
for ib in bonds:
s = self._S[ib]
if len(s.shape) == 1:
res.append(entropy(s**2, n))
else:
if for_matrix_S:
# explicitly calculate Schmidt values by diagonalizing (s^dagger s)
s = npc.eigvalsh(npc.tensordot(s.conj(), s, axes=[0, 0]))
res.append(entropy(s, n))
else:
raise ValueError("entropy with non-diagonal schmidt values")
return np.array(res)
[docs] def entanglement_entropy_segment(self, segment=[0], first_site=None, n=1):
r"""Calculate entanglement entropy for general geometry of the bipartition.
This function is similar as :meth:`entanglement_entropy`,
but for more general geometry of the region `A` to be a segment of a *few* sites.
This is acchieved by explicitly calculating the reduced density matrix of `A`
and thus works only for small segments.
Parameters
----------
segment : list of int
Given a first site `i`, the region ``A_i`` is defined to be ``[i+j for j in segment]``.
first_site : ``None`` | (iterable of) int
Calculate the entropy for segments starting at these sites.
``None`` defaults to ``range(L-segment[-1])`` for finite
or `range(L)` for infinite boundary conditions.
n : int | float
Selects which entropy to calculate;
`n=1` (default) is the ususal von-Neumann entanglement entropy,
otherwise the `n`-th Renyi entropy.
Returns
-------
entropies : 1D ndarray
``entropies[i]`` contains the entropy for the the region ``A_i`` defined above.
"""
# Side-Remark: there is a trick to calculate the entanglement for large regions `A_i`
# of consecutive sites (in our notation, ``segment = range(La)``)
# To get the entanglement entropy, diagonalize:
# --theta---
# | | |
# --theta*--
# Diagonalization is O(chi^6), compared to O(d^{3*La})
segment = np.sort(segment)
if first_site is None:
if self.finite:
first_site = range(0, self.L - segment[-1])
else:
first_site = range(self.L)
comb_legs = [
self._get_p_labels(len(segment), False),
self._get_p_labels(len(segment), True)
]
res = []
for i0 in first_site:
rho = self.get_rho_segment(segment + i0)
rho = rho.combine_legs(comb_legs, qconj=[+1, -1])
p = npc.eigvalsh(rho)
res.append(entropy(p, n))
return np.array(res)
[docs] def entanglement_spectrum(self, by_charge=False):
r"""return entanglement energy spectrum.
Parameters
----------
by_charge : bool
Wheter we should sort the spectrum on each bond by the possible charges.
Returns
-------
ent_spectrum : list
For each (non-trivial) bond the entanglement spectrum.
If `by_charge` is ``False``, return (for each bond) a sorted 1D ndarray
with the convention :math:`S_i^2 = e^{-\xi_i}`, where :math:`S_i` labels a Schmidt
value and :math:`\xi_i` labels the entanglement 'energy' in the returned spectrum.
If `by_charge` is True, return a a list of tuples ``(charge, sub_spectrum)``
for each possible charge on that bond.
"""
if by_charge:
res = []
for i in range(self.L + 1)[self.nontrivial_bonds]:
ss = -2. * np.log(self._S[i])
if i < self.L:
leg = self._B[i].get_leg('vL')
else: # i == L: segment b.c.
leg = self._B[i - 1].get_leg('vR').conj()
spectrum = [(leg.get_charge(qi), np.sort(ss[leg.get_slice(qi)]))
for qi in range(leg.block_number)]
res.append(spectrum)
return res
else:
return [np.sort(-2. * np.log(ss)) for ss in self._S[self.nontrivial_bonds]]
[docs] def get_rho_segment(self, segment):
"""Return reduced density matrix for a segment.
Note that the dimension of rho_A scales exponentially in the length of the segment.
Parameters
----------
segment : iterable of int
Sites for which the reduced density matrix is to be calculated.
Assumed to be sorted.
Returns
-------
rho : :class:`~tenpy.linalg.np_conserved.Array`
Reduced density matrix of the segment sites.
Labels ``'p0', 'p1', ..., 'pk', 'p0*', 'p1*', ..., 'pk*'`` with ``k=len(segment)``.
"""
if len(segment) > 12:
warnings.warn("{0:d} sites in the segment, that's much!".format(len(segment)),
stacklevel=2)
if len(segment) > 20:
raise ValueError("too large segment; this is exponentially expensive!")
segment = np.sort(segment)
if np.all(segment[1:] == segment[:-1] + 1): # consecutive
theta = self.get_theta(segment[0], segment[-1] - segment[0] + 1)
rho = npc.tensordot(theta, theta.conj(), axes=(['vL', 'vR'], ['vL*', 'vR*']))
return rho
rho = self.get_theta(segment[0], 1)
rho = npc.tensordot(rho, rho.conj(), axes=('vL', 'vL*'))
k = 1
contract_axes = (['vR*'] + self._p_label, ['vL*'] + self._get_p_label('*'))
for i in range(segment[0] + 1, segment[-1]):
B = self.get_B(i)
if i == segment[k]:
B = self._replace_p_label(B, str(k))
k += 1
rho = npc.tensordot(rho, B, axes=('vR', 'vL'))
rho = npc.tensordot(rho, B.conj(), axes=('vR*', 'vL*'))
else:
rho = npc.tensordot(rho, B, axes=('vR', 'vL'))
rho = npc.tensordot(rho, B.conj(), axes=contract_axes)
B = self._replace_p_label(self.get_B(segment[-1]), str(k))
rho = npc.tensordot(rho, B, axes=('vR', 'vL'))
rho = npc.tensordot(rho, B.conj(), axes=(['vR*', 'vR'], ['vL*', 'vR*']))
return rho
[docs] def probability_per_charge(self, bond=0):
"""Return probabilites of charge value on the left of a given bond.
For example for particle number conservation, define
:math:`N_b = sum_{i<b} n_i` for a given bond `b`.
This function returns the possible values of `N_b` as rows of `charge_values`,
and for each row the probabilty that this combination occurs in the given state.
Parameters
----------
bond : int
The bond to be considered. The returned charges are summed on the left of this bond.
Returns
-------
charge_values : 2D array
Columns correspond to the different charges in `self.chinfo`.
Rows are the different charge fluctuations at this bond
probabilities : 1D array
For each row of `charge_values` the probablity for these values of charge fluctuations.
"""
if self.bc == 'segment' and bond == self.L:
S = self.get_SR(self.L - 1)**2
leg = self.get_B(self.L - 1, form=None).get_leg('vR').conj()
else: # usually the case
S = self.get_SL(bond)**2
leg = self.get_B(bond, form=None).get_leg('vL')
assert leg.qconj == +1
if not leg.is_blocked():
raise ValueError("leg not blocked: can have duplicate entries in charge values")
ps = []
for qi in range(leg.block_number):
sl = leg.get_slice(qi)
ps.append(np.sum(S[sl]))
ps = np.array(ps)
if abs(np.sum(ps) - 1.) > 1.e-10:
warnings.warn("Probability_per_charge: Sum of probabilites not 1. Canonical form?",
stacklevel=2)
return leg.charges.copy(), ps
[docs] def average_charge(self, bond=0):
r"""Return the average charge for the block on the left of a given bond.
For example for particle number conservation, define
:math:`N_b = sum_{i<b} n_i` for a given bond `b`.
Then this function returns :math:`<\psi| N_b |\psi>`.
Parameters
----------
bond : int
The bond to be considered.
The returned charges are summed over the sites left of `bond`.
Returns
-------
average_charge : 1D array
For each type of charge in :attr:`chinfo`
the average value when summing the charge values over sites left of the given bond.
"""
charges, ps = self.probability_per_charge(bond)
return np.sum(ps[:, np.newaxis] * charges, axis=0)
[docs] def charge_variance(self, bond=0):
r"""Return the charge variance on the left of a given bond.
For example for particle number conservation, define
:math:`N_b = sum_{i<b} n_i` for a given bond `b`.
Then this function returns :math:`<\psi| N_b^2 |\psi> - (<\psi| N_b |\psi>)^2`.
Parameters
----------
bond : int
The bond to be considered.
The returned charges are summed over the sites left of `bond`.
Returns
-------
average_charge : 1D array
For each type of charge in :attr:`chinfo`
the variance of of the charge values left of the given bond.
"""
charges_mean = self.average_charge(bond)
charges, ps = self.probability_per_charge(bond)
return np.sum(ps[:, np.newaxis] * (charges - charges_mean[:, np.newaxis])**2, axis=0)
[docs] def mutinf_two_site(self, max_range=None, n=1):
"""Calculate the two-site mutual information :math:`I(i:j)`.
Calculates :math:`I(i:j) = S(i) + S(j) - S(i,j)`,
where :math:`S(i)` is the single site entropy on site :math:`i`
and :math:`S(i,j)` the two-site entropy on sites :math:`i,j`.
Parameters
----------
max_range : int
Maximal distance ``|i-j|`` for which the mutual information should be calculated.
``None`` defaults to `L-1`.
n : float
Selects the entropy to use, see :func:`~tenpy.tools.math.entropy`.
Returns
-------
coords : 2D array
Coordinates for the mutinf array.
mutinf : 1D array
``mutinf[k]`` is the mutual information :math:`I(i:j)` between the
sites ``i, j = coords[k]``.
"""
# Basically the code of get_rho_segment and entanglement_entropy,
# but optimized to run in O(L*max_range)
if max_range is None:
max_range = self.L
S_i = self.entanglement_entropy_segment(n=n) # single-site entropy
legs_ij = self._get_p_labels(2, False), self._get_p_labels(2, True)
# = (['p0', 'p1'], ['p0*', 'p1*'])
contr_legs = (
['vR*'] + self._get_p_label('1'), # ['vL', 'p1']
['vL*'] + self._get_p_label('1*')) # ['vL*', 'p1*']
mutinf = []
coord = []
for i in range(self.L):
rho = self.get_theta(i, 1)
rho = npc.tensordot(rho, rho.conj(), axes=('vL', 'vL*'))
jmax = i + max_range + 1
if self.finite:
jmax = min(jmax, self.L)
for j in range(i + 1, jmax):
B = self._replace_p_label(self.get_B(j, form='B'), '1') # 'vL', 'vR', 'p1'
rho = npc.tensordot(rho, B, axes=['vR', 'vL'])
rho_ij = npc.tensordot(rho, B.conj(), axes=(['vR*', 'vR'], ['vL*', 'vR*']))
rho_ij = rho_ij.combine_legs(legs_ij, qconj=[+1, -1])
S_ij = entropy(npc.eigvalsh(rho_ij), n)
mutinf.append(S_i[i] + S_i[j % self.L] - S_ij)
coord.append((i, j))
if j + 1 < jmax:
rho = npc.tensordot(rho, B.conj(), axes=contr_legs)
return np.array(coord), np.array(mutinf)
[docs] def overlap(self, other, charge_sector=0, ignore_form=False, **kwargs):
"""Compute overlap ``<self|other>``.
Parameters
----------
other : :class:`MPS`
An MPS with the same physical sites.
charge_sector : None | charges | ``0``
Selects the charge sector in which the dominant eigenvector of the TransferMatrix is.
``None`` stands for *all* sectors, ``0`` stands for the zero-charge sector.
Defaults to ``0``, i.e., *assumes* the dominant eigenvector is in charge sector 0.
ignore_form : bool
If ``False`` (default), take into account the canonical form :attr:`form` at each site.
If ``True``, we ignore the canonical form (i.e., whether the MPS is in left, right,
mixed or no canonical form) and just contract all the :attr:`_B` as they are.
(This can give different results!)
**kwargs :
Further keyword arguments given to :meth:`TransferMatrix.eigenvectors`;
only used for infinite boundary conditions.
Returns
-------
overlap : dtype.type
The contraction ``<self|other> * self.norm * other.norm``
(i.e., taking into account the :attr:`norm` of both MPS).
For an infinite MPS, ``<self|other>`` is the overlap per unit cell, i.e.,
the largest eigenvalue of the TransferMatrix.
"""
if self.finite:
if ignore_form:
# Use TransferMatrix with option to ignore the form
TM = TransferMatrix(self, other, charge_sector=charge_sector, form=None)
res = TM.matvec(TM.initial_guess(1.)) # apply transfer matrix to identity
return npc.trace(res, 0, 1) * self.norm * other.norm
else:
env = MPSEnvironment(self, other)
return env.full_contraction(0)
else: # infinite
form = None if ignore_form else 'B'
TM = TransferMatrix(self, other, charge_sector=charge_sector, form=form)
ov, _ = TM.eigenvectors(**kwargs)
return ov[0] * self.norm * other.norm
[docs] def expectation_value(self, ops, sites=None, axes=None):
"""Expectation value ``<psi|ops|psi>/<psi|psi>`` of (n-site) operator(s).
Given the MPS in canonical form, it calculates n-site expectation values.
For example the contraction for a two-site (`n` = 2) operator on site `i` would look like::
| .--S--B[i]--B[i+1]--.
| | | | |
| | |-----| |
| | | op | |
| | |-----| |
| | | | |
| .--S--B*[i]-B*[i+1]-.
Parameters
----------
ops : (list of) { :class:`~tenpy.linalg.np_conserved.Array` | str }
The operators, for wich the expectation value should be taken,
All operators should all have the same number of legs (namely `2 n`).
If less than `self.L` operators are given, we repeat them periodically.
Strings (like ``'Id', 'Sz'``) are translated into single-site operators defined by
:attr:`sites`.
sites : None | list of int
List of site indices. Expectation values are evaluated there.
If ``None`` (default), the entire chain is taken (clipping for finite b.c.)
axes : None | (list of str, list of str)
Two lists of each `n` leg labels giving the physical legs of the operator used for
contraction. The first `n` legs are contracted with conjugated `B`,
the second `n` legs with the non-conjugated `B`.
``None`` defaults to ``(['p'], ['p*'])`` for single site operators (`n` = 1), or
``(['p0', 'p1', ... 'p{n-1}'], ['p0*', 'p1*', .... 'p{n-1}*'])`` for `n` > 1.
Returns
-------
exp_vals : 1D ndarray
Expectation values, ``exp_vals[i] = <psi|ops[i]|psi>``, where ``ops[i]`` acts on
site(s) ``j, j+1, ..., j+{n-1}`` with ``j=sites[i]``.
Examples
--------
One site examples (n=1):
>>> psi.expectation_value('Sz')
[Sz0, Sz1, ..., Sz{L-1}]
>>> psi.expectation_value(['Sz', 'Sx'])
[Sz0, Sx1, Sz2, Sx3, ... ]
>>> psi.expectation_value('Sz', sites=[0, 3, 4])
[Sz0, Sz3, Sz4]
Two site example (n=2), assuming homogeneous sites:
>>> SzSx = npc.outer(psi.sites[0].Sz.replace_labels(['p', 'p*'], ['p0', 'p0*']),
psi.sites[1].Sx.replace_labels(['p', 'p*'], ['p1', 'p1*']))
>>> psi.expectation_value(SzSx)
[Sz0Sx1, Sz1Sx2, Sz2Sx3, ... ] # with len L-1 for finite bc, or L for infinite
Example measuring <psi|SzSx|psi2> on each second site, for inhomogeneous sites:
>>> SzSx_list = [npc.outer(psi.sites[i].Sz.replace_labels(['p', 'p*'], ['p0', 'p0*']),
psi.sites[i+1].Sx.replace_labels(['p', 'p*'], ['p1', 'p1*']))
for i in range(0, psi.L-1, 2)]
>>> psi.expectation_value(SzSx_list, range(0, psi.L-1, 2))
[Sz0Sx1, Sz2Sx3, Sz4Sx5, ...]
"""
ops, sites, n, (op_ax_p, op_ax_pstar) = self._expectation_value_args(ops, sites, axes)
ax_p = ['p' + str(k) for k in range(n)]
ax_pstar = ['p' + str(k) + '*' for k in range(n)]
E = []
for i in sites:
op = self.get_op(ops, i)
op = op.replace_labels(op_ax_p + op_ax_pstar, ax_p + ax_pstar)
theta = self.get_theta(i, n)
C = npc.tensordot(op, theta, axes=[ax_pstar, ax_p]) # C has same labels as theta
E.append(npc.inner(theta, C, axes='labels', do_conj=True))
return np.real_if_close(np.array(E))
[docs] def expectation_value_term(self, term, autoJW=True):
r"""Expectation value ``<psi|op_{i0}op_{i1}...op_{iN}|psi>/<psi|psi>``.
Calculates the expectation value of a tensor product of single-site operators
acting on different sites `i0`, `i1`, ... (not necessarily next to each other).
In other words, evaluate the expectation value of a term ``op0_i0 op1_i1 op2_i2 ...``.
For example the contraction of three one-site operators on sites `i0`,
`i1=i0+1`, `i2=i0+3` would look like::
| .--S--B[i0]---B[i0+1]--B[i0+2]--B[i0+3]--.
| | | | | | |
| | op1 op2 | op3 |
| | | | | | |
| .--S--B*[i0]--B*[i0+1]-B*[i0+2]-B*[i0+3]-.
Parameters
----------
term : list of (str, int)
List of tuples ``op, i`` where `i` is the MPS index of the site the operator
named `op` acts on.
The order inside `term` determines the order in which they act
(in the mathematical convention: the last operator in `term` is right-most,
so it acts first on a Ket).
autoJW : bool
If True (default), automatically insert Jordan Wigner strings for Fermions as needed.
Returns
-------
exp_val : float/complex
The expectation value of the tensorproduct of the given onsite operators,
``<psi|op_i0 op_i1 ... op_iN |psi>/<psi|psi>``,
where ``|psi>`` is the represented MPS.
See also
--------
correlation_function : efficient way to evaluate many correlation functions.
Examples
--------
>>> a = psi.expectation_value_term([('Sx', 2), ('Sz', 4)])
>>> b = psi.expectation_value_term([('Sz', 4), ('Sx', 2)])
>>> c = psi.expectation_value_multi_sites(['Sz', 'Id', 'Sz'], i0=2)
>>> assert a == b == c
"""
# strategy: translate term into a list "ops" to be used for `expectation_value_multi_sites`
term = list(term)
i_min = min([t[1] for t in term])
i_max = max([t[1] for t in term])
ops = [None] * (i_max - i_min + 1)
count_JW = 0
for op, i in term:
j = i - i_min # index in ops
if ops[j] is not None:
ops[j] = ops[j] + " " + op
else:
ops[j] = op
if autoJW and self.sites[self._to_valid_index(i)].op_needs_JW(op):
count_JW += 1
for k in range(j):
if ops[k] is not None:
ops[k] = ops[k] + ' JW'
if ops[k].endswith(' JW JW'):
ops[k] = ops[k][:-len(' JW JW')]
else:
ops[k] = 'JW'
for i in range(len(ops)):
if ops[i] is None:
ops[i] = 'Id'
if count_JW % 2 == 1:
raise ValueError("Odd number of operators which need a Jordan Wigner string")
return self.expectation_value_multi_sites(ops, i_min)
[docs] def expectation_value_multi_sites(self, operators, i0):
r"""Expectation value ``<psi|op0_{i0}op1_{i0+1}...opN_{i0+N}|psi>/<psi|psi>``.
Calculates the expectation value of a tensor product of single-site operators
acting on different sites next to each other.
In other words, evaluate the expectation value of a term
``op0_i0 op1_{i0+1} op2_{i0+2} ...``.
Parameters
----------
operators : List of { :class:`~tenpy.linalg.np_conserved.Array` | str }
List of one-site operators. This method calculates the
expectation value of the n-sites operator given by their tensor
product.
i0 : int
The left most index on which an operator acts, i.e.,
``operators[i]`` acts on site ``i + i0``.
Returns
-------
exp_val : float/complex
The expectation value of the tensorproduct of the given onsite operators,
``<psi|operators[0]_{i0} operators[1]_{i0+1} ... |psi>/<psi|psi>``,
where ``|psi>`` is the represented MPS.
"""
op = operators[0]
if (isinstance(op, str)):
op = self.sites[self._to_valid_index(i0)].get_op(op)
theta = self.get_B(i0, 'Th')
C = npc.tensordot(op, theta, axes=['p*', 'p'])
axes = [['vL*'] + self._get_p_label('*'), ['vL'] + self._p_label]
C = npc.tensordot(theta.conj(), C, axes=axes)
axes[1][0] = 'vR*'
for j in range(1, len(operators)):
op = operators[j] # the operator
i = i0 + j # the site it acts on
B = self.get_B(i, form='B')
C = npc.tensordot(C, B, axes=['vR', 'vL'])
if op != 'Id':
if (isinstance(op, str)):
op = self.sites[self._to_valid_index(i)].get_op(op)
C = npc.tensordot(op, C, axes=['p*', 'p'])
C = npc.tensordot(B.conj(), C, axes=axes)
exp_val = npc.trace(C, 'vR*', 'vR')
return np.real_if_close(exp_val)
[docs] def expectation_value_terms_sum(self, term_list, prefactors=None):
"""Calculate expectation values for a bunch of terms and sum them up.
This is equivalent to the following expression::
sum([self.expectation_value_term(term)*strength for term, strength in term_list])
However, for effiency, the term_list is converted to an MPO and the expectation value
of the MPO is evaluated.
.. note ::
Due to the way MPO expectation values are evaluated for infinite systems,
it works only if all terms in the `term_list` start within the MPS unit cell.
.. deprecated:: 0.4.0
`prefactor` will be removed in version 1.0.0.
Instead, directly give just ``TermList(term_list, prefactors)`` as argument.
Parameters
----------
term_list : :class:`~tenpy.networks.terms.TermList`
The terms and prefactors (`strength`) to be summed up.
prefactors :
Instead of specifying a :class:`~tenpy.networks.terms.TermList`,
one can also specify the term_list and strength separately.
This is deprecated.
Returns
-------
terms_sum : list of (complex) float
Equivalent to the expression
``sum([self.expectation_value_term(term)*strength for term, strength in term_list])``.
_mpo :
Intermediate results: the generated MPO.
For a finite MPS, ``terms_sum = _mpo.expectation_value(self)``, for an infinite MPS
``terms_sum = _mpo.expectation_value(self) * self.L``
See also
--------
expectation_value_term : evaluates a single `term`.
tenpy.networks.mpo.MPO.expectation_value : expectation value density of an MPO.
"""
from . import mpo, terms
if prefactors is not None:
warnings.warn(
"Deprecated argument prefactors: replace arguments with "
"``TermList(term_list, prefactors)``.", FutureWarning, 2)
term_list = terms.TermList(term_list, prefactors)
L = self.L
if not self.finite:
for term in term_list.terms:
if not 0 <= min([i for _, i in term]) < L:
raise ValueError("term doesn't start in MPS unit cell: " + repr(term))
# conversion
ot, ct = term_list.to_OnsiteTerms_CouplingTerms(self.sites)
bc = 'finite' if self.finite else 'infinite'
mpo_graph = mpo.MPOGraph.from_terms(ot, ct, self.sites, bc)
mpo_ = mpo_graph.build_MPO()
terms_sum = mpo_.expectation_value(self, max_range=ct.max_range())
if not self.finite:
terms_sum = terms_sum * self.L
return terms_sum, mpo_
[docs] def correlation_function(self,
ops1,
ops2,
sites1=None,
sites2=None,
opstr=None,
str_on_first=True,
hermitian=False):
r"""Correlation function ``<psi|op1_i op2_j|psi>/<psi|psi>`` of single site operators.
Given the MPS in canonical form, it calculates 2-site correlation functions.
For examples the contraction for a two-site operator on site `i` would look like::
| .--S--B[i]--B[i+1]--...--B[j]---.
| | | | | |
| | | | op2 |
| | op1 | | |
| | | | | |
| .--S--B*[i]-B*[i+1]-...--B*[j]--.
Onsite terms are taken in the order ``<psi | op1 op2 | psi>``.
If `opstr` is given and ``str_on_first=True``, it calculates::
| for i < j for i > j
|
| .--S--B[i]---B[i+1]--...- B[j]---. .--S--B[j]---B[j+1]--...- B[i]---.
| | | | | | | | | | |
| | opstr opstr op2 | | op2 | | |
| | | | | | | | | | |
| | op1 | | | | opstr opstr op1 |
| | | | | | | | | | |
| .--S--B*[i]--B*[i+1]-...- B*[j]--. .--S--B*[j]--B*[j+1]-...- B*[i]--.
For ``i==j``, no `opstr` is included.
For ``str_on_first=False``, the `opstr` on site ``min(i, j)`` is always left out.
Strings (like ``'Id', 'Sz'``) in the arguments are translated into single-site
operators defined by the :class:`~tenpy.networks.site.Site` on which they act.
Each operator should have the two legs ``'p', 'p*'``.
Parameters
----------
ops1 : (list of) { :class:`~tenpy.linalg.np_conserved.Array` | str }
First operator of the correlation function (acting after ops2).
``ops1[x]`` acts on site ``sites1[x]``.
If less than ``len(sites1)`` operators are given, we repeat them periodically.
ops2 : (list of) { :class:`~tenpy.linalg.np_conserved.Array` | str }
Second operator of the correlation function (acting before ops1).
``ops2[y]`` acts on site ``sites2[y]``.
If less than ``len(sites2)`` operators are given, we repeat them periodically.
sites1 : None | int | list of int
List of site indices; a single `int` is translated to ``range(0, sites1)``.
``None`` defaults to all sites ``range(0, L)``.
Is sorted before use, i.e. the order is ignored.
sites2 : None | int | list of int
List of site indices; a single `int` is translated to ``range(0, sites2)``.
``None`` defaults to all sites ``range(0, L)``.
Is sorted before use, i.e. the order is ignored.
opstr : None | (list of) { :class:`~tenpy.linalg.np_conserved.Array` | str }
Ignored by default (``None``).
Operator(s) to be inserted between ``ops1`` and ``ops2``.
If less than :attr:`L` operators are given, we repeat them periodically.
If given as a list, ``opstr[r]`` is inserted at site `r` (independent of `sites1` and
`sites2`).
str_on_first : bool
Whether the `opstr` is included on the site ``min(i, j)``.
Note the order, which is chosen that way to handle fermionic Jordan-Wigner strings
correctly. (In other words: choose ``str_on_first=True`` for fermions!)
hermitian : bool
Optimization flag: if ``sites1 == sites2`` and ``Ops1[i]^\dagger == Ops2[i]``
(which is not checked explicitly!), the resulting ``C[x, y]`` will be hermitian.
We can use that to avoid calculations, so ``hermitian=True`` will run faster.
Returns
-------
C : 2D ndarray
The correlation function ``C[x, y] = <psi|ops1[i] ops2[j]|psi>``,
where ``ops1[i]`` acts on site ``i=sites1[x]`` and ``ops2[j]`` on site ``j=sites2[y]``.
If `opstr` is given, it gives (for ``str_on_first=True``):
- For ``i < j``: ``C[x, y] = <psi|ops1[i] prod_{i <= r < j} opstr[r] ops2[j]|psi>``.
- For ``i > j``: ``C[x, y] = <psi|prod_{j <= r < i} opstr[r] ops1[i] ops2[j]|psi>``.
- For ``i = j``: ``C[x, y] = <psi|ops1[i] ops2[j]|psi>``.
The condition ``<= r`` is replaced by a strict ``< r``, if ``str_on_first=False``.
"""
ops1, ops2, sites1, sites2, opstr = self._correlation_function_args(
ops1, ops2, sites1, sites2, opstr)
if hermitian and sites1 != sites2:
warnings.warn("MPS correlation function can't use the hermitian flag", stacklevel=2)
hermitian = False
C = np.empty((len(sites1), len(sites2)), dtype=np.complex)
for x, i in enumerate(sites1):
# j > i
j_gtr = sites2[sites2 > i]
if len(j_gtr) > 0:
C_gtr = self._corr_up_diag(ops1, ops2, i, j_gtr, opstr, str_on_first, True)
C[x, (sites2 > i)] = C_gtr
if hermitian:
C[x + 1:, x] = np.conj(C_gtr)
# j == i
j_eq = sites2[sites2 == i]
if len(j_eq) > 0:
# on-site correlation function
op12 = npc.tensordot(self.get_op(ops1, i), self.get_op(ops2, i), axes=['p*', 'p'])
C[x, (sites2 == i)] = self.expectation_value(op12, i, [['p'], ['p*']])
if not hermitian:
# j < i
for y, j in enumerate(sites2):
i_gtr = sites1[sites1 > j]
if len(i_gtr) > 0:
C[(sites1 > j), y] = self._corr_up_diag(ops2, ops1, j, i_gtr, opstr,
str_on_first, False)
# exchange ops1 and ops2 : they commute on different sites,
# but we apply opstr after op1 (using the last argument = False)
return np.real_if_close(C)
[docs] def norm_test(self):
"""Check that self is in canonical form.
Returns
-------
norm_error: array, shape (L, 2)
For each site the norm error to the left and right.
The error ``norm_error[i, 0]`` is defined as the norm-difference between
the following networks::
| --theta[i]---. --s[i]--.
| | | vs |
| --theta*[i]--. --s[i]--.
Similarly, ``norm_errror[i, 1]`` is the norm-difference of::
| .--theta[i]--- .--s[i+1]--
| | | vs |
| .--theta*[i]-- .--s[i+1]--
"""
err = np.empty((self.L, 2), dtype=np.float)
lbl_R = (self._get_p_label('0') + ['vR'], self._get_p_label('0*') + ['vR*'])
lbl_L = (['vL'] + self._get_p_label('0'), ['vL*'] + self._get_p_label('0*'))
for i in range(self.L):
th = self.get_theta(i, 1)
rho_L = npc.tensordot(th, th.conj(), axes=lbl_R)
S = self.get_SL(i)
if isinstance(S, npc.Array): # during DMRG with mixer, S may be a 2D npc.Array
if S.rank != 2:
raise ValueError("Expect 2D npc.Array or 1D numpy ndarray")
rho_L2 = npc.tensordot(S, S.conj(), axes=['vR', 'vR*'])
else:
rho_L2 = npc.diag(S**2, rho_L.get_leg('vL'), dtype=rho_L.dtype)
err[i, 0] = npc.norm(rho_L - rho_L2)
rho_R = npc.tensordot(th, th.conj(), axes=lbl_L)
S = self.get_SR(i)
if isinstance(S, npc.Array):
if S.rank != 2:
raise ValueError("Expect 2D npc.Array or 1D numpy ndarray")
rho_R2 = npc.tensordot(S, S.conj(), axes=['vL', 'vL*'])
else:
rho_R2 = npc.diag(S**2, rho_R.get_leg('vR'), dtype=rho_L.dtype)
err[i, 1] = npc.norm(rho_R - rho_R2)
return err
[docs] def correlation_length(self, target=1, tol_ev0=1.e-8, charge_sector=0):
r"""Calculate the correlation length by diagonalizing the transfer matrix.
Assumes that `self` is in canonical form.
Works only for infinite MPS, where the transfer matrix is a useful concept.
Assuming a single-site unit cell, any correlation function splits into
:math:`C(A_i, B_j) = A'_i T^{j-i-1} B'_j`
with some parts left and right and the :math:`j-i-1`-th power of the transfer matrix in
between. The largest eigenvalue is 1 (if self is properly normalized)
and gives the dominant contribution of
:math:`A'_i E_1 * 1^{j-i-1} * E_1^T B'_j = <A> <B>`,
and the second largest one gives a contribution :math:`\propto \lambda_2^{j-i-1}`.
Thus :math:`\lambda_2 = \exp(-\frac{1}{\xi})`.
More general for a `L`-site unit cell we get :math:`\lambda_2 = \exp(-\frac{L}{\xi})`,
where the `xi` is given in units of 1 lattice spacing in the MPS.
.. warning ::
For a higher-dimensional lattice (which the MPS class doesn't know about),
the correct unit is the lattice spacing in x-direction, and the correct formula is
:math:`\lambda_2 = \exp(-\frac{L_x}{\xi})`,
where `L_x` is the number of lattice spacings in the infinite direction within the
MPS unit cell, e.g. the number of "rings" of a cylinder in the MPS unit cell.
To get to these units, divide the returned `xi` by the number of sites within a "ring",
for a lattice given in :attr:`~tenpy.networks.lattice.N_sites_per_ring`.
Parameters
----------
target : int
We look for the `target` + 1 largest eigenvalues.
tol_ev0 : float
Print warning if largest eigenvalue deviates from 1 by more than `tol_ev0`.
charge_sector : None | charges | ``0``
Selects the charge sector in which the dominant eigenvector of the TransferMatrix is.
``None`` stands for *all* sectors, ``0`` stands for the zero-charge sector.
Defaults to ``0``, i.e., *assumes* the dominant eigenvector is in charge sector 0.
Returns
-------
xi : float | 1D array
If `target`=1, return just the correlation length,
otherwise an array of the `target` largest correlation lengths.
It is measured in units of a single lattice spacing in the MPS language,
see the warning above.
"""
assert (not self.finite)
T = TransferMatrix(self, self, charge_sector=charge_sector, form='B')
num = max(target + 1, self._transfermatrix_keep)
E, _ = T.eigenvectors(num, which='LM')
E = E[np.argsort(-np.abs(E))] # sort descending by magnitude
if charge_sector is not None and charge_sector != 0:
# need also dominant eigenvector: include 0 charge sector to results
del T
T = TransferMatrix(self, self, charge_sector=0, form='B')
E0, _ = T.eigenvectors(num, which='LM')
assert abs(E0[0]) > abs(E[0]), "dominant eigenvector in zero charge sector?"
E = np.array([E0[0]] + list(E))
if abs(E[0] - 1.) > tol_ev0:
warnings.warn(
"Correlation length: largest eigenvalue not one. "
"Not in canonical form/normalized?",
stacklevel=2)
if len(E) < 2:
return 0. # only a single eigenvector: zero correlation length
if target == 1:
return -1. / np.log(abs(E[1] / E[0])) * self.L
return -1. / np.log(np.abs(E[1:target + 1] / E[0])) * self.L
[docs] def add(self, other, alpha, beta, cutoff=1.e-15):
"""Return an MPS which represents ``alpha|self> + beta |others>``.
Works only for 'finite', 'segment' boundary conditions.
For 'segment' boundary conditions, the virtual legs on the very left/right are
assumed to correspond to each other (i.e. self and other have the same state outside of
the considered segment).
Takes into account :attr:`norm`.
Parameters
----------
other : :class:`MPS`
Another MPS of the same length to be added with self.
alpha, beta : complex float
Prefactors for self and other. We calculate
``alpha * |self> + beta * |other>``
cutoff : float | None
Cutoff of singular values used in the SVDs.
Returns
-------
sum : :class:`MPS`
An MPS representing ``alpha|self> + beta |other>``.
Has same total charge as `self`.
U_L, V_R : :class:`~tenpy.linalg.np_conserved.Array`
Only returned for ``'segment'`` boundary conditions.
The unitaries defining the new left and right Schmidt states in terms of the old ones,
with legs ``'vL', 'vR'``.
"""
L = self.L
assert (other.L == L and L >= 2) # (if you need this, generalize this function...)
assert self.finite
self._gauge_compatible_vL_vR(other)
legs = ['vL', 'vR'] + self._p_label
# alpha and beta appear only on the first site
alpha = alpha * self.norm
beta = beta * other.norm
Bs = [
npc.grid_concat(
[[alpha * self.get_B(0).transpose(legs), beta * other.get_B(0).transpose(legs)]],
axes=[0, 1])
]
for i in range(1, L - 1):
B1 = self.get_B(i).transpose(legs)
B2 = other.get_B(i).transpose(legs)
grid = [[B1, npc.zeros([B1.get_leg('vL'), B2.get_leg('vR')] + B1.legs[2:])],
[npc.zeros([B2.get_leg('vL'), B1.get_leg('vR')] + B1.legs[2:]), B2]]
Bs.append(npc.grid_concat(grid, [0, 1]))
Bs.append(
npc.grid_concat(
[[self.get_B(L - 1).transpose(legs)], [other.get_B(L - 1).transpose(legs)]],
axes=[0, 1]))
Ss = [np.ones(1)] + [np.ones(B.shape[1]) for B in Bs]
psi = self.__class__(self.sites, Bs, Ss, 'finite', form=None) # new class instance
# bring to canonical form, calculate Ss
if self.bc == 'segment':
U_L, V_R = psi.canonical_form_finite(renormalize=False, cutoff=cutoff)
return psi, U_L, V_R
else:
psi.canonical_form_finite(renormalize=False, cutoff=cutoff)
return psi
[docs] def apply_local_op(self, i, op, unitary=None, renormalize=False, cutoff=1.e-13):
"""Apply a local (one or multi-site) operator to `self`.
Note that this destroys the canonical form if the local operator is non-unitary.
Therefore, this function calls :meth:`canonical_form` if necessary.
Parameters
----------
i : int
(Left-most) index of the site(s) on which the operator should act.
op : str | npc.Array
A physical operator acting on site `i`, with legs ``'p', 'p*'`` for a single-site
operator or with legs ``['p0', 'p1', ...], ['p0*', 'p1*', ...]`` for an operator
acting on `n`>=2 sites.
Strings (like ``'Id', 'Sz'``) are translated into single-site operators defined by
:attr:`sites`.
unitary : None | bool
Whether `op` is unitary, i.e., whether the canonical form is preserved (``True``)
or whether we should call :meth:`canonical_form` (``False``).
``None`` checks whether ``norm(op dagger(op) - identity)`` is smaller than `cutoff`.
renormalize : bool
Whether the final state should keep track of the norm (False, default) or be
renormalized to have norm 1 (True).
cutoff : float
Cutoff for singular values if `op` acts on more than one site (see :meth:`from_full`).
(And used as cutoff for a unspecified `unitary`.)
"""
i = self._to_valid_index(i)
if isinstance(op, str):
op = self.sites[i].get_op(op)
n = op.rank // 2 # same as int(rank/2)
if n == 1:
pstar, p = 'p*', 'p'
else:
p = self._get_p_labels(n, False)
pstar = self._get_p_labels(n, True)
if unitary is None:
op_op_dagger = npc.tensordot(op, op.conj(), axes=[pstar, p])
if n > 1:
op_op_dagger = op_op_dagger.combine_legs([p, pstar], qconj=[+1, -1])
unitary = npc.norm(op_op_dagger - npc.eye_like(op_op_dagger)) < cutoff
if n == 1:
opB = npc.tensordot(op, self._B[i], axes=['p*', 'p'])
self.set_B(i, opB, self.form[i])
else:
th = self.get_theta(i, n)
th = npc.tensordot(op, th, axes=[pstar, p])
# use MPS.from_full to split the sites
split_th = self.from_full(self.sites[i:i + n], th, None, cutoff, False, 'segment',
(self.get_SL(i), self.get_SR(i + n - 1)))
for j in range(n):
self.set_B(i + j, split_th._B[j], split_th.form[j])
for j in range(n - 1):
self.set_SR(i + j, split_th._S[j + 1])
if not unitary:
self.canonical_form(renormalize)
[docs] def swap_sites(self, i, swap_op='auto', trunc_par={}):
"""Swap the two neighboring sites `i` and `i+1` (inplace).
Exchange two neighboring sites: form theta, 'swap' the physical legs and split
with an svd. While the 'swap' is just a transposition/relabeling for bosons, one needs to
be careful about the sign for fermions.
Parameters
----------
i : int
Swap the two sites at positions `i` and `i+1`.
swap_op : ``None`` | ``'auto'`` | :class:`~tenpy.linalg.np_conserved.Array`
The operator used to swap the phyiscal legs of the two-site wave function `theta`.
For ``None``, just transpose/relabel the legs, for ``'auto'`` also take care of
fermionic signs. Alternative give an npc :class:`~tenpy.linalg.np_conserved.Array`
which represents the full operator used for the swap.
Should have legs ``['p0', 'p1', 'p0*', 'p1*']`` whith ``'p0', 'p1*'`` contractible.
trunc_par : dict
Parameters for truncation, see :func:`~tenpy.algorithms.truncation.truncate`.
`chi_max` defaults to ``max(self.chi)``.
Returns
-------
trunc_err : :class:`~tenpy.algorithms.truncation.TruncationError`
The error of the represented state introduced by the truncation after the swap.
"""
trunc_par.setdefault('chi_max', max(self.chi))
siteL, siteR = self.sites[self._to_valid_index(i)], self.sites[self._to_valid_index(i + 1)]
if swap_op == 'auto':
# get sign for Fermions.
# If we write the wave function as
# psi = sum_{ [n_i]} psi_[n_i] prod_i (c^dagger_i)^{n_i} |vac>
# we see that switching i <-> i+1 the phase to be introduced is by commuting
# (c^dagger_i)^{n_i} with (c^dagger_{i+1})^{n_{i+1}}
# This gives a sign (-1)^{n_i * n_{i+1}}.
# site.JW_exponent is the `n_i` in the above equations, for each physical index.
sign = siteL.JW_exponent[:, np.newaxis] * siteR.JW_exponent[np.newaxis, :]
if np.any(sign):
dL, dR = siteL.dim, siteR.dim
sign = np.real_if_close(np.exp(1.j * np.pi * sign.reshape([dL * dR])))
swap_op = np.diag(sign).reshape([dL, dR, dL, dR])
legs = [siteL.leg, siteR.leg, siteL.leg.conj(), siteR.leg.conj()]
swap_op = npc.Array.from_ndarray(swap_op, legs)
swap_op.iset_leg_labels(['p1', 'p0', 'p0*', 'p1*'])
else: # no sign necessary
swap_op = None # continue with transposition as for Bosons
theta = self.get_theta(i, n=2)
C = self.get_theta(i, n=2, formL=0.) # inversion free, see also TEBDEngine.update_bond()
if swap_op is None:
# just replace the labels, effectively this is a transposition.
theta.ireplace_labels(['p0', 'p1'], ['p1', 'p0'])
C.ireplace_labels(['p0', 'p1'], ['p1', 'p0'])
elif isinstance(swap_op, npc.Array):
theta = npc.tensordot(swap_op, theta, axes=[['p0*', 'p1*'], ['p0', 'p1']])
C = npc.tensordot(swap_op, C, axes=(['p0*', 'p1*'], ['p0', 'p1']))
else:
raise ValueError("Invalid swap_op: got " + repr(swap_op))
theta = theta.combine_legs([('vL', 'p0'), ('vR', 'p1')], qconj=[+1, -1])
U, S, V, err, renormalize = svd_theta(theta, trunc_par, inner_labels=['vR', 'vL'])
B_R = V.split_legs(1).ireplace_label('p1', 'p')
B_L = npc.tensordot(C.combine_legs(('vR', 'p1'), pipes=theta.legs[1]),
V.conj(),
axes=['(vR.p1)', '(vR*.p1*)'])
B_L.ireplace_labels(['vL*', 'p0'], ['vR', 'p'])
B_L /= renormalize # re-normalize to <psi|psi> = 1
self.set_SR(i, S)
self.set_B(i, B_L, 'B')
self.set_B(i + 1, B_R, 'B')
self.sites[self._to_valid_index(i)] = siteR # swap 'sites' as well
self.sites[self._to_valid_index(i + 1)] = siteL
return err
[docs] def permute_sites(self, perm, swap_op='auto', trunc_par={}, verbose=0):
"""Applies the permutation perm to the state (inplace).
Parameters
----------
perm : ndarray[ndim=1, int]
The applied permutation, such that ``psi.permute_sites(perm)[i] = psi[perm[i]]``
(where ``[i]`` indicates the `i`-th site).
swap_op : ``None`` | ``'auto'`` | :class:`~tenpy.linalg.np_conserved.Array`
The operator used to swap the phyiscal legs of a two-site wave function `theta`,
see :meth:`swap_sites`.
trunc_par : dict
Parameters for truncation, see :func:`~tenpy.algorithms.truncation.truncate`.
`chi_max` defaults to ``max(self.chi)``.
verbose : float
Level of verbosity, print status messages if verbose > 0.
Returns
-------
trunc_err : :class:`~tenpy.algorithms.truncation.TruncationError`
The error of the represented state introduced by the truncation after the swaps.
"""
# In order to keep sites close together, we always scan from the left,
# keeping everything up to `i` in strictly ascending order.
# => more or less an 'insertion' sort algorithm.
# Works nicely for permutations like [1,2,3,0,6,7,8,5] (swapping the 0 and 5 around).
# For [ 2 3 4 5 6 7 0 1], it splits 0 and 1 apart (first swapping the 0 down, then the 1)
trunc_par.setdefault('chi_max', max(self.chi))
trunc_err = TruncationError()
num_swaps = 0
i = 0
while i < self.L - 1:
if perm[i] > perm[i + 1]:
if verbose > 1:
print(i, ": chi = ", self._S[i + 1].shape[0], end='')
trunc = self.swap_sites(i, swap_op, trunc_par)
if verbose > 1:
print("->", self._S[i + 1].shape[0], ". eps = ", trunc.eps)
num_swaps += 1
x, y = perm[i], perm[i + 1]
perm[i + 1], perm[i] = x, y
# restart from very left; but we know it's already sorted up to i-1
if i > 0:
i -= 1
trunc_err += trunc
else:
i += 1
if verbose > 0:
print("Total swaps in permute_sites:", num_swaps, repr(trunc_err))
return trunc_err
[docs] def compute_K(self, perm, swap_op='auto', trunc_par=None, canonicalize=1.e-6, verbose=0):
r"""Compute the momentum quantum numbers of the entanglement spectrum for 2D states.
Works for an infinite MPS living on a cylinder, infinitely long in `x` direction and with
periodic boundary conditions in `y` directions.
If the state is invariant under 'rotations' around the cylinder axis, one can find the
momentum quantum numbers of it. (The rotation is nothing more than a translation in `y`.)
This function permutes some sites (on a copy of `self`) to enact the rotation, and then
finds the dominant eigenvector of the mixed transfer matrix to get the quantum numbers,
along the lines of [PollmannTurner2012]_, see also (the appendix and Fig. 11 in the arXiv
version of) [CincioVidal2013]_.
Parameters
----------
perm : 1D ndarray | :class:`~tenpy.models.lattice.Lattice`
Permuation to be applied to the physical indices, see :meth:`permute_sites`.
If a lattice is given, we use it to read out the lattice structure and shift
each site by one lattice-vector in y-direction (assuming periodic boundary conditions).
(If you have a :class:`~tenpy.models.model.CouplingModel`,
give its `lat` attribute for this argument)
swap_op : ``None`` | ``'auto'`` | :class:`~tenpy.linalg.np_conserved.Array`
The operator used to swap the phyiscal legs of a two-site wave function `theta`,
see :meth:`swap_sites`.
trunc_par : dict
Parameters for truncation, see :func:`~tenpy.algorithms.truncation.truncate`.
If not set, `chi_max` defaults to ``max(self.chi)``.
canonicalize : float
Check that `self` is in canonical form; call :meth:`canonical_form`
if :meth:`norm_test` yields ``np.linalg.norm(self.norm_test()) > canonicalize``.
verbose : float
Level of verbosity, print status messages if verbose > 0.
Returns
-------
U : :class:`~tenpy.linalg.np_conserved.Array`
Unitary representation of the applied permutation on left Schmidt states.
W : ndarray
1D array of the form ``S**2 exp(i K)``, where `S` are the Schmidt values
on the left bond. You can use :func:`np.abs` and :func:`np.angle` to extract the
Schmidt values `S` and momenta `K` from `W`.
q : :class:`~tenpy.linalg.charges.LegCharge`
LegCharge corresponding to `W`.
ov : complex
The eigenvalue of the mixed transfer matrix `<psi|T|psi>` per :attr:`L` sites.
An absolute value different smaller than 1 indicates that the state is not invariant
under the permutation or that the truncation error `trunc_err` was too large!
trunc_err : :class:`~tenpy.algorithms.truncation.TruncationError`
The error of the represented state introduced by the truncation after swaps when
performing the truncation.
"""
from ..models.lattice import Lattice # dynamical import to avoid import loops
if self.finite:
raise ValueError("Works only for infinite b.c.")
if trunc_par is None:
trunc_par = {}
trunc_par.setdefault('chi_max', max(self.chi))
trunc_par.setdefault('verbose', verbose)
if isinstance(perm, Lattice):
lat = perm
assert lat.dim >= 2 # ensure that the lattice is at least 2D
assert lat.N_sites == self.L
shifted_lat_order = lat.order.copy()
shifted_lat_order[:, 1] = np.mod(shifted_lat_order[:, 1] + 1, lat.Ls[1])
perm = lat.lat2mps_idx(shifted_lat_order)
if verbose > 1:
print("permutation: ", perm)
# preliminary: check canonical form
self.convert_form('B')
norm_err = np.linalg.norm(self.norm_test())
if norm_err > canonicalize:
warnings.warn("self.norm_test() = {0!s} ==> canonicalize".format(self.norm_test()))
self.canonical_form()
# get copy of self
psi_t = self.copy()
# apply permutation
perm = np.asarray(perm)
trunc_err = psi_t.permute_sites(perm, swap_op, trunc_par, verbose / 10.)
# re-check canonical form
norm_err = np.linalg.norm(psi_t.norm_test())
if norm_err > canonicalize:
warnings.warn("psi_t.norm_test() = {0!s} ==> canonicalize".format(psi_t.norm_test()))
psi_t.convert_form('B')
TM = TransferMatrix(self, psi_t, transpose=True, charge_sector=0)
# Find left dominant eigenvector of this mixed transfer matrix.
# Because we are in B form and get the left eigenvector,
# the resulting vector should be sUs up to a scaling.
ov, sUs = TM.eigenvectors(num_ev=self._transfermatrix_keep)
if verbose > 0:
print("compute_K: overlap ", ov[0], ", |o| = 1. -", 1. - np.abs(ov[0]))
# (should be 1 if state is invariant under translations)
print("compute_K: truncation error ", trunc_err.eps)
sUs = sUs[0].split_legs(0)
_, sUs_blocked = sUs.as_completely_blocked()
W = npc.eigvals(sUs_blocked, sort='m>')
# W = s^2 exp(i K ) up to overall scaling
# Strip S's from U
inv_S = 1. / self.get_SL(0)
U = sUs.scale_axis(inv_S, 0).iscale_axis(inv_S, 1)
# U should be unitary - scale it
U *= (np.sqrt(U.shape[0]) / npc.norm(U))
return U, W / np.sum(np.abs(W)), sUs_blocked.legs[0], ov[0], trunc_err
def __str__(self):
"""Some status information about the MPS."""
res = ["MPS, L={L:d}, bc={bc!r}.".format(L=self.L, bc=self.bc)]
res.append("chi: " + str(self.chi))
if self.L > 10:
res.append("first two sites: " + repr(self.sites[0]) + " " + repr(self.sites[1]))
res.append("first two forms:" + " ".join([repr(f) for f in self.form[:2]]))
else:
res.append("sites: " + " ".join([repr(s) for s in self.sites]))
res.append("forms: " + " ".join([repr(f) for f in self.form]))
return "\n".join(res)
def _to_valid_index(self, i):
"""Make sure `i` is a valid index (depending on `self.bc`)."""
if not self.finite:
return i % self.L
if i < 0:
i += self.L
if i >= self.L or i < 0:
raise ValueError("i = {0:d} out of bounds for finite MPS".format(i))
return i
def _parse_form(self, form):
"""Parse `form` = (list of) {tuple | key of _valid_forms} to list of tuples"""
if isinstance(form, tuple):
return [form] * self.L
form = to_iterable(form)
if len(form) == 1:
form = [form[0]] * self.L
if len(form) != self.L:
raise ValueError("Wrong len of `form`: " + repr(form))
return [self._to_valid_form(f) for f in form]
def _to_valid_form(self, form):
"""Parse `form` = {tuple | key of _valid_forms} to a tuple"""
if isinstance(form, tuple):
return form
return self._valid_forms[form]
def _scale_axis_B(self, B, S, form_diff, axis_B, cutoff):
"""Scale an axis of B with S to bring it in desired form.
If S is just 1D (as usual, e.g. during TEBD), this function just performs
``B.scale_axis(S**form_diff, axis_B)``.
However, during the DMRG with mixer, S might acutally be a 2D matrix.
For ``form_diff = -1``, we need to calculate the inverse of S, more precisely the
(Moore-Penrose) pseudo inverse, see :func:`~tenpy.linalg.np_conserved.pinv`.
The cutoff is only used in that case.
Returns scaled B.
"""
if form_diff == 0:
return B # nothing to do
if not isinstance(S, npc.Array):
# the usual case: S is a 1D array with singular values
if form_diff != 1.:
S = S**form_diff
return B.scale_axis(S, axis_B)
else:
# e.g. during DMRG with a DensityMatrixMixer
if S.rank != 2:
raise ValueError("Expect 2D npc.Array or 1D numpy ndarray")
if form_diff == -1:
S = npc.pinv(S, cutoff)
elif form_diff != 1.:
raise ValueError("Can't scale/tensordot a 2D `S` for non-integer `form_diff`")
# Hack: mpo.MPOEnvironment.full_contraction uses ``axis_B == 'vL*'``
if axis_B == 'vL' or axis_B == 'vL*':
B = npc.tensordot(S, B, axes=[1, axis_B]).replace_label(0, axis_B)
elif axis_B == 'vR' or axis_B == 'vR*':
B = npc.tensordot(B, S, axes=[axis_B, 0]).replace_label(-1, axis_B)
else:
raise ValueError("This should never happen: unexpected leg for scaling with S")
return B
def _replace_p_label(self, A, s):
"""Return npc Array `A` with replaced label, ``'p' -> 'p'+s``.
This is done for each of the 'physical labels' in :attr:`_p_label`. With a clever use of
this function, the re-implementation of various functions (like get_theta) in derived
classes with multiple legs per site can be avoided.
"""
return A.replace_label('p', 'p' + s)
# return A.replace_labels(self._p_label, self._get_p_label(s))
def _get_p_label(self, s):
"""return self._p_label with additional string `s`."""
return ['p' + s]
# return [lbl + s for lbl in self._p_label]
def _get_p_labels(self, ks, star=False):
"""Join ``self._get_p_label(str(k)) for k in range(ks)`` to a single list."""
if star:
return ['p' + str(k) + '*' for k in range(ks)]
# return [lbl + str(k) + '*' for k in range(ks) for lbl in self._p_label]
else:
return ['p' + str(k) for k in range(ks)]
# return [lbl + str(k) for k in range(ks) for lbl in self._p_label]
def _expectation_value_args(self, ops, sites, axes):
"""parse the arguments of self.expectation_value()"""
ops = npc.to_iterable_arrays(ops)
if any(isinstance(op, str) for op in ops):
n = 1
else:
s = 0 if sites is None else to_iterable(sites)[0]
n = ops[s % len(ops)].rank // 2 # same as int(rank/2)
L = self.L
if sites is None:
if self.finite:
sites = range(L - (n - 1))
else:
sites = range(L)
sites = to_iterable(sites)
if axes is None:
if n == 1:
axes = (['p'], ['p*'])
else:
axes = (self._get_p_labels(n), self._get_p_labels(n, True))
# check number of axes
ax_p, ax_pstar = axes
if len(ax_p) != n or len(ax_pstar) != n:
raise ValueError("Len of axes does not match to n-site operator with n=" + str(n))
return ops, sites, n, axes
def _correlation_function_args(self, ops1, ops2, sites1, sites2, opstr):
"""get default arguments of self.correlation_function()"""
if sites1 is None:
sites1 = range(0, self.L)
elif isinstance(sites1, int):
sites1 = range(0, sites1)
if sites2 is None:
sites2 = range(0, self.L)
elif isinstance(sites2, int):
sites2 = range(0, sites2)
ops1 = npc.to_iterable_arrays(ops1)
ops2 = npc.to_iterable_arrays(ops2)
opstr = npc.to_iterable_arrays(opstr)
sites1 = np.sort(sites1)
sites2 = np.sort(sites2)
return ops1, ops2, sites1, sites2, opstr
def _corr_up_diag(self, ops1, ops2, i, j_gtr, opstr, str_on_first, apply_opstr_first):
"""correlation function above the diagonal: for fixed i and all j in j_gtr, j > i."""
op1 = self.get_op(ops1, i)
opstr1 = self.get_op(opstr, i)
if opstr1 is not None:
axes = ['p*', 'p'] if apply_opstr_first else ['p', 'p*']
op1 = npc.tensordot(op1, opstr1, axes=axes)
theta = self.get_theta(i, n=1)
C = npc.tensordot(op1, theta, axes=['p*', 'p0'])
C = npc.tensordot(theta.conj(), C, axes=[['p0*', 'vL*'], ['p', 'vL']])
# C has legs 'vR*', 'vR'
js = list(j_gtr[::-1]) # stack of j, sorted *descending*
res = []
for r in range(i + 1, js[0] + 1): # js[0] is the maximum
B = self.get_B(r, form='B')
C = npc.tensordot(C, B, axes=['vR', 'vL'])
if r == js[-1]:
Cij = npc.tensordot(self.get_op(ops2, r), C, axes=['p*', 'p'])
Cij = npc.inner(B.conj(), Cij, axes=[['vL*', 'p*', 'vR*'], ['vR*', 'p', 'vR']])
res.append(Cij)
js.pop()
if len(js) > 0:
op = self.get_op(opstr, r)
if op is not None:
C = npc.tensordot(op, C, axes=['p*', 'p'])
C = npc.tensordot(B.conj(), C, axes=[['vL*', 'p*'], ['vR*', 'p']])
return res
def _canonical_form_dominant_gram_matrix(self, bond0, transpose, tol_xi, guess=None):
"""Find dominant eigenvector of the transfer matrix starting between sites (bond0-1,bond0).
Find right (transpose=False) or left (transpose=True) eigenvector of the transfermatrix.
"""
TM = TransferMatrix(self, self, bond0, transpose=transpose, charge_sector=0)
if guess is None:
diag = self.get_SL(bond0)**2 if transpose else 1.
guess = TM.initial_guess(diag)
guess = guess.combine_legs([0, 1], pipes=TM.pipe)
eta, V = TM.eigenvectors(self._transfermatrix_keep, v0=guess, which='LM')
self._transfermatrix_keep = len(eta)
if len(eta) > 1:
if np.abs(eta[0]) > np.abs(eta[1]):
xi = -self.L / np.log(np.abs((eta[1] / eta[0])))
else:
xi = np.inf
if xi > tol_xi:
raise ValueError("Degenerate spectrum of TransferMatrix "
"(corr length xi={xi:.3e})".format(xi=xi))
eta, G = eta[0], V[0]
G = G.split_legs()
# note: the dominant eigenvector should be hermitian and positive
# removes phase (arbitrary for eigenvectors!) and normalize
# right eigenvectors should have trace chi, left ones trace 1
# (since we expect something close to eye(chi) for right and diag(S**2) for left G)
norm = 1. if transpose else G.shape[0]
G *= norm / npc.trace(G, 0, 1)
if self.dtype.kind != 'c': # psi is real -> G should be real
eta = np.abs(eta)
G.iunary_blockwise(np.real)
return eta, G # G has legs vL, vL* or vR, vR*
def _canonical_form_correct_right(self,
i1,
Gr,
return_Gl_guess=False,
eps=2. * np.finfo(np.double).eps):
"""Given the right gram matrix Gr, updated the bond (i0, i1), where i0 = i1 - 1.
Diagonalise Gr = X^H Wr X and update
``B[i0] -> B[i0] X^H / norm``,
``B[i1] -> X B[i1] * norm``, where norm = sqrt(chi/sum(Wr)) == sqrt(chi/tr(Gr))
If `Wr` has (almost) zero entries, reduce the bond dimension at the given bond.
Then ``Gr -> Wr``.
Return Wr normalized to ``sum(Wr) = chi``.
If `return_Gl_guess`, return also ``Gl_guess = X S[i1]**2 X^H``.
"""
Gr.itranspose(['vL', 'vL*'])
W, XH = npc.eigh(Gr) # -> XH has legs vL vL* = vL vR
if np.sign(W[np.argmax(np.abs(W))]) == -1: # fix sign
W = -W # should actually never happen: we initially normalize tr(Gr) = chi > 0
# discard small values on order of machine precision
proj = (W > eps)
if np.count_nonzero(proj) < len(W):
# project into non-degenerate subspace, reducing the bond dimensions!
warnings.warn("canonical_form_infinite: project to smaller bond dimension",
stacklevel=3)
XH.iproject(proj, axes=1)
W = W[proj]
norm = len(W) / np.sum(W)
W *= norm
norm = np.sqrt(norm) # (norm doesn't change eigenvalue of TM)
Kl = XH.iset_leg_labels(['vL', 'vR']) * (1. / norm)
Kr = XH.transpose().iconj().iset_leg_labels(['vL', 'vR']) * norm
i0 = i1 - 1
self.set_B(i0, npc.tensordot(self.get_B(i0), Kl, axes=['vR', 'vL']))
self.set_B(i1, npc.tensordot(Kr, self.get_B(i1), axes=['vR', 'vL']))
if return_Gl_guess:
guess = npc.tensordot(Kr.scale_axis(self.get_SL(i1)**2, 1), Kl, axes=['vR', 'vL'])
return W, guess.iset_leg_labels(['vR*', 'vR'])
return W
def _canonical_form_correct_left(self, i1, Gl, Wr, eps=2. * np.finfo(np.double).eps):
"""Bring into canonical form on bond (i0, i1) where i0= i1 - 1.
Given the left Gram matrix Gl (with legs 'vR*', 'vR')
and right diag(Wr), compute and diagonalize the density matrix
``rho = sqrt(Wr) Gl sqrt(Wr) -> Y^H S^2 Y`` (Y acting on ket, Y^H on bra).
Then we can update
B[i0] -> B[i0] sqrt(Wr) Y^H
B[i1] -> Y 1/sqrt(Wr) B[i1]
Thus the new dominant left eigenvector is
Gl -> Y sqrt(Wr) Gl sqrt(Wr) Y^H = S^2
and dominant right eigenvector is
diag(Wr) -> Y 1/sqrt(Wr) diag(Wr) 1/sqrt(W) Y^H = diag(1)
i.e., we brought the bond to canonical form and `S` is the Schmidt spectrum.
"""
sqrt_Wr = np.sqrt(Wr)
Gl.itranspose(['vR*', 'vR'])
rhor = Gl.scale_axis(sqrt_Wr, 0).iscale_axis(sqrt_Wr, 1)
S2, YH = npc.eigh(rhor, sort='>') # YH has legs 'vR*', 'vR'
S2 /= np.sum(S2) # equivalent to normalizing tr(rhor)=1
s_norm = 1.
# discard small values on order of machine precision
proj = (S2 > eps)
if np.count_nonzero(proj) < len(S2):
# project into non-degenerate subspace, reducing the bond dimensions!
warnings.warn("canonical_form_infinite: project to smaller bond dimension",
stacklevel=2)
YH.iproject(proj, axes=1)
S2 = S2[proj]
s_norm = np.sqrt(np.sum(S2))
S = np.sqrt(S2) / s_norm
self.set_SL(i1, S)
Yl = YH.scale_axis(sqrt_Wr / s_norm, 0).iset_leg_labels(['vL', 'vR'])
Yr = YH.transpose().iconj().scale_axis(1. / sqrt_Wr, 1).iset_leg_labels(['vL', 'vR'])
i0 = i1 - 1
self.set_B(i0, npc.tensordot(self.get_B(i0), Yl, axes=['vR', 'vL']))
self.set_B(i1, npc.tensordot(Yr, self.get_B(i1), axes=['vR', 'vL']))
Gl = npc.tensordot(Gl, Yl, axes=['vR', 'vL'])
Gl = npc.tensordot(Yl.conj(), Gl, axes=['vL*', 'vR*']) # labels 'vR*', 'vR'
Gl /= npc.trace(Gl)
# Gl is diag(S**2) up to numerical errors...
return Gl, np.ones(Yr.legs[0].ind_len, np.float)
def _gauge_compatible_vL_vR(self, other):
"""If necessary, gauge total charge of `other` to match the vL, vR legs of self."""
if self.chinfo.qnumber == 0:
return
from tenpy.tools import optimization
need_gauge = False
with optimization.temporary_level(optimization.OptimizationFlag.default):
try:
other._B[0].get_leg('vL').test_equal(self._B[0].get_leg('vL'))
other._B[-1].get_leg('vR').test_equal(self._B[-1].get_leg('vR'))
except ValueError:
need_gauge = True
if need_gauge:
other.gauge_total_charge(None, self._B[0].get_leg('vL'), self._B[-1].get_leg('vR'))
if any(self.get_total_charge() != other.get_total_charge()):
raise ValueError("self and other have different total charges!")
[docs]class MPSEnvironment:
"""Stores partial contractions of :math:`<bra|Op|ket>` for local operators `Op`.
The network for a contraction :math:`<bra|Op|ket>` of a local operator `Op`, say exemplary
at sites `i, i+1` looks like::
| .-----M[0]--- ... --M[1]---M[2]--- ... ->--.
| | | | | |
| | | |------| |
| LP[0] | | Op | RP[-1]
| | | |------| |
| | | | | |
| .-----N[0]*-- ... --N[1]*--N[2]*-- ... -<--.
Of course, we can also calculate the overlap `<bra|ket>` by using the special case ``Op = Id``.
We use the following label convention (where arrows indicate `qconj`)::
| .-->- vR vL ->-.
| | |
| LP RP
| | |
| .--<- vR* vL* -<-.
To avoid recalculations of the whole network e.g. in the DMRG sweeps,
we store the contractions up to some site index in this class.
For ``bc='finite','segment'``, the very left and right part ``LP[0]`` and
``RP[-1]`` are trivial and don't change,
but for ``bc='infinite'`` they are might be updated
(by inserting another unit cell to the left/right).
The MPS `bra` and `ket` have to be in canonical form.
All the environments are constructed without the singular values on the open bond.
In other words, we contract left-canonical `A` to the left parts `LP`
and right-canonical `B` to the right parts `RP`.
Thus, the special case ``ket=bra`` should yield identity matrices for `LP` and `RP`.
Parameters
----------
bra : :class:`~tenpy.networks.mps.MPS`
The MPS to project on. Should be given in usual 'ket' form;
we call `conj()` on the matrices directly.
Stored in place, without making copies.
If necessary to match charges, we call :meth:`~tenpy.networks.mps.MPS.gauge_total_charge`.
ket : :class:`~tenpy.networks.mpo.MPO` | None
The MPS on which the local operator acts.
Stored in place, without making copies.
If ``None``, use `bra`.
init_LP : ``None`` | :class:`~tenpy.linalg.np_conserved.Array`
Initial very left part ``LP``. If ``None``, build trivial one with :meth:`init_LP`.
init_RP : ``None`` | :class:`~tenpy.linalg.np_conserved.Array`
Initial very right part ``RP``. If ``None``, build trivial one with :meth:`init_RP`.
age_LP : int
The number of physical sites involved into the contraction yielding `firstLP`.
age_RP : int
The number of physical sites involved into the contraction yielding `lastRP`.
Attributes
----------
L : int
Number of physical sites involved into the Environment, i.e. the least common multiple
of ``bra.L`` and ``ket.L``.
bra, ket : :class:`~tenpy.networks.mps.MPS`
The two MPS for the contraction.
dtype : type
The data type.
_finite : bool
Whether the boundary conditions of the MPS are finite.
_LP : list of {``None`` | :class:`~tenpy.linalg.np_conserved.Array`}
Left parts of the environment, len `L`.
``LP[i]`` contains the contraction strictly left of site `i`
(or ``None``, if we don't have it calculated).
_RP : list of {``None`` | :class:`~tenpy.linalg.np_conserved.Array`}
Right parts of the environment, len `L`.
``RP[i]`` contains the contraction strictly right of site `i`
(or ``None``, if we don't have it calculated).
_LP_age : list of int | ``None``
Used for book-keeping, how large the DMRG system grew:
``_LP_age[i]`` stores the number of physical sites invovled into the contraction
network which yields ``self._LP[i]``.
_RP_age : list of int | ``None``
Used for book-keeping, how large the DMRG system grew:
``_RP_age[i]`` stores the number of physical sites invovled into the contraction
network which yields ``self._RP[i]``.
"""
def __init__(self, bra, ket, init_LP=None, init_RP=None, age_LP=0, age_RP=0):
if ket is None:
ket = bra
if ket is not bra:
ket._gauge_compatible_vL_vR(bra) # ensure matching charges
self.bra = bra
self.ket = ket
self.dtype = np.find_common_type([bra.dtype, ket.dtype], [])
self.L = L = lcm(bra.L, ket.L)
self._finite = bra.finite
self._LP = [None] * L
self._RP = [None] * L
self._LP_age = [None] * L
self._RP_age = [None] * L
self._finite = self.ket.finite # just for _to_valid_index
if init_LP is None:
init_LP = self.init_LP(0)
self.set_LP(0, init_LP, age=age_LP)
if init_RP is None:
init_RP = self.init_RP(L - 1)
self.set_RP(L - 1, init_RP, age=age_RP)
self.test_sanity()
[docs] def test_sanity(self):
"""Sanity check, raises ValueErrors, if something is wrong."""
assert (self.bra.finite == self.ket.finite == self._finite)
# check that the network is contractable
for i in range(self.L):
b_s = self.bra.sites[i % self.bra.L]
k_s = self.ket.sites[i % self.ket.L]
b_s.leg.test_equal(k_s.leg)
assert any([LP is not None for LP in self._LP])
assert any([RP is not None for RP in self._RP])
[docs] def init_LP(self, i):
"""Build initial left part ``LP``.
Parameters
----------
i : int
Build ``LP`` left of site `i`.
Returns
-------
init_LP : :class:`~tenpy.linalg.np_conserved.Array`
Identity contractible with the `vL` leg of ``ket.get_B(i)``, labels ``'vR*', 'vR'``.
"""
leg_ket = self.ket.get_B(i, None).get_leg('vL')
leg_bra = self.bra.get_B(i, None).get_leg('vL')
leg_ket.test_equal(leg_bra)
init_LP = npc.diag(1., leg_ket, dtype=self.dtype)
init_LP.iset_leg_labels(['vR*', 'vR'])
return init_LP
[docs] def init_RP(self, i):
"""Build initial right part ``RP`` for an MPS/MPOEnvironment.
Parameters
----------
i : int
Build ``RP`` right of site `i`.
Returns
-------
init_RP : :class:`~tenpy.linalg.np_conserved.Array`
Identity contractible with the `vR` leg of ``ket.get_B(i)``, labels ``'vL*', 'vL'``.
"""
leg_ket = self.ket.get_B(i, None).get_leg('vR')
leg_bra = self.bra.get_B(i, None).get_leg('vR')
leg_ket.test_equal(leg_bra)
init_RP = npc.diag(1., leg_ket, dtype=self.dtype)
init_RP.iset_leg_labels(['vL*', 'vL'])
return init_RP
[docs] def get_LP(self, i, store=True):
"""Calculate LP at given site from nearest available one (including `i`).
The returned ``LP_i`` corresponds to the following contraction,
where the M's and the N's are in the 'A' form::
| .-------M[0]--- ... --M[i-1]--->- 'vR'
| | | |
| LP[0] | |
| | | |
| .-------N[0]*-- ... --N[i-1]*--<- 'vR*'
Parameters
----------
i : int
The returned `LP` will contain the contraction *strictly* left of site `i`.
store : bool
Wheter to store the calculated `LP` in `self` (``True``) or discard them (``False``).
Returns
-------
LP_i : :class:`~tenpy.linalg.np_conserved.Array`
Contraction of everything left of site `i`,
with labels ``'vR*', 'vR'`` for `bra`, `ket`.
"""
# find nearest available LP to the left.
for i0 in range(i, i - self.L, -1):
LP = self._LP[self._to_valid_index(i0)]
if LP is not None:
break
# (for finite, LP[0] should always be set, so we should abort at latest with i0=0)
else: # no break called
raise ValueError("No left part in the system???")
age_i0 = self.get_LP_age(i0)
for j in range(i0, i):
LP = self._contract_LP(j, LP)
if store:
self.set_LP(j + 1, LP, age=age_i0 + j - i0 + 1)
return LP
[docs] def get_RP(self, i, store=True):
"""Calculate RP at given site from nearest available one (including `i`).
The returned ``RP_i`` corresponds to the following contraction,
where the M's and the N's are in the 'B' form::
| 'vL' ->---M[i+1]-- ... --M[L-1]----.
| | | |
| | | RP[-1]
| | | |
| 'vL*' -<---N[i+1]*- ... --N[L-1]*---.
Parameters
----------
i : int
The returned `RP` will contain the contraction *strictly* right of site `i`.
store : bool
Wheter to store the calculated `RP` in `self` (``True``) or discard them (``False``).
Returns
-------
RP_i : :class:`~tenpy.linalg.np_conserved.Array`
Contraction of everything left of site `i`,
with labels ``'vL*', 'vL'`` for `bra`, `ket`.
"""
# find nearest available RP to the right.
for i0 in range(i, i + self.L):
RP = self._RP[self._to_valid_index(i0)]
if RP is not None:
break
# (for finite, RP[-1] should always be set, so we should abort at latest with i0=L-1)
age_i0 = self.get_RP_age(i0)
for j in range(i0, i, -1):
RP = self._contract_RP(j, RP)
if store:
self.set_RP(j - 1, RP, age=age_i0 + i0 - j + 1)
return RP
[docs] def get_LP_age(self, i):
"""Return number of physical sites in the contractions of get_LP(i).
Might be ``None``.
"""
return self._LP_age[self._to_valid_index(i)]
[docs] def get_RP_age(self, i):
"""Return number of physical sites in the contractions of get_RP(i).
Might be ``None``.
"""
return self._RP_age[self._to_valid_index(i)]
[docs] def set_LP(self, i, LP, age):
"""Store part to the left of site `i`."""
i = self._to_valid_index(i)
self._LP[i] = LP
self._LP_age[i] = age
[docs] def set_RP(self, i, RP, age):
"""Store part to the right of site `i`."""
i = self._to_valid_index(i)
self._RP[i] = RP
self._RP_age[i] = age
[docs] def del_LP(self, i):
"""Delete stored part strictly to the left of site `i`."""
i = self._to_valid_index(i)
self._LP[i] = None
self._LP_age[i] = None
[docs] def del_RP(self, i):
"""Delete storde part scrictly to the right of site `i`."""
i = self._to_valid_index(i)
self._RP[i] = None
self._RP_age[i] = None
[docs] def full_contraction(self, i0):
"""Calculate the overlap by a full contraction of the network.
The full contraction of the environments gives the overlap ``<bra|ket>``,
taking into account :attr:`MPS.norm` of both `bra` and `ket`.
For this purpose, this function contracts
``get_LP(i0+1, store=False)`` and ``get_RP(i0, store=False)`` with appropriate singular
values in between.
Parameters
----------
i0 : int
Site index.
"""
if self.ket.finite and i0 + 1 == self.L:
# special case to handle `_to_valid_index` correctly:
# get_LP(L) is not valid for finite b.c, so we use need to calculate it explicitly.
LP = self.get_LP(i0, store=False)
LP = self._contract_LP(i0, LP)
else:
LP = self.get_LP(i0 + 1, store=False)
# multiply with `S`: a bit of a hack: use 'private' MPS._scale_axis_B
S_bra = self.bra.get_SR(i0).conj()
LP = self.bra._scale_axis_B(LP, S_bra, form_diff=1., axis_B='vR*', cutoff=0.)
# cutoff is not used for form_diff = 1
S_ket = self.ket.get_SR(i0)
LP = self.bra._scale_axis_B(LP, S_ket, form_diff=1., axis_B='vR', cutoff=0.)
RP = self.get_RP(i0, store=False)
contr = npc.inner(LP, RP, axes=[['vR*', 'vR'], ['vL*', 'vL']], do_conj=False)
return contr * self.bra.norm * self.ket.norm
[docs] def expectation_value(self, ops, sites=None, axes=None):
"""Expectation value ``<bra|ops|ket>`` of (n-site) operator(s).
Calculates n-site expectation values of operators sandwiched between bra and ket.
For examples the contraction for a two-site operator on site `i` would look like::
| .--S--B[i]--B[i+1]--.
| | | | |
| | |-----| |
| LP[i] | op | RP[i+1]
| | |-----| |
| | | | |
| .--S--B*[i]-B*[i+1]-.
Here, the `B` are taken from `ket`, the `B*` from `bra`.
The call structure is the same as for :meth:`MPS.expectation_value`.
.. warning :
In contrast to :meth:`MPS.expectation_value`, this funciton does not normalize,
thus it also takes into account :attr:`MPS.norm` of both `bra` and `ket`.
Parameters
----------
ops : (list of) { :class:`~tenpy.linalg.np_conserved.Array` | str }
The operators, for wich the expectation value should be taken,
All operators should all have the same number of legs (namely `2 n`).
If less than ``len(sites)`` operators are given, we repeat them periodically.
Strings (like ``'Id', 'Sz'``) are translated into single-site operators defined by
:attr:`sites`.
sites : list
List of site indices. Expectation values are evaluated there.
If ``None`` (default), the entire chain is taken (clipping for finite b.c.)
axes : None | (list of str, list of str)
Two lists of each `n` leg labels giving the physical legs of the operator used for
contraction. The first `n` legs are contracted with conjugated `B`,
the second `n` legs with the non-conjugated `B`.
``None`` defaults to ``(['p'], ['p*'])`` for single site (n=1), or
``(['p0', 'p1', ... 'p{n-1}'], ['p0*', 'p1*', .... 'p{n-1}*'])`` for `n` > 1.
Returns
-------
exp_vals : 1D ndarray
Expectation values, ``exp_vals[i] = <bra|ops[i]|ket>``, where ``ops[i]`` acts on
site(s) ``j, j+1, ..., j+{n-1}`` with ``j=sites[i]``.
Examples
--------
One site examples (n=1):
>>> env.expectation_value('Sz')
[Sz0, Sz1, ..., Sz{L-1}]
>>> env.expectation_value(['Sz', 'Sx'])
[Sz0, Sx1, Sz2, Sx3, ... ]
>>> env.expectation_value('Sz', sites=[0, 3, 4])
[Sz0, Sz3, Sz4]
Two site example (n=2), assuming homogeneous sites:
>>> SzSx = npc.outer(psi.sites[0].Sz.replace_labels(['p', 'p*'], ['p0', 'p0*']),
psi.sites[1].Sx.replace_labels(['p', 'p*'], ['p1', 'p1*']))
>>> env.expectation_value(SzSx)
[Sz0Sx1, Sz1Sx2, Sz2Sx3, ... ] # with len L-1 for finite bc, or L for infinite
Example measuring <bra|SzSx|ket> on each second site, for inhomogeneous sites:
>>> SzSx_list = [npc.outer(psi.sites[i].Sz.replace_labels(['p', 'p*'], ['p0', 'p0*']),
psi.sites[i+1].Sx.replace_labels(['p', 'p*'], ['p1', 'p1*']))
for i in range(0, psi.L-1, 2)]
>>> env.expectation_value(SzSx_list, range(0, psi.L-1, 2))
[Sz0Sx1, Sz2Sx3, Sz4Sx5, ...]
"""
ops, sites, n, (op_ax_p, op_ax_pstar) = self.ket._expectation_value_args(ops, sites, axes)
ax_p = ['p' + str(k) for k in range(n)]
ax_pstar = ['p' + str(k) + '*' for k in range(n)]
E = []
for i in sites:
LP = self.get_LP(i, store=True)
RP = self.get_RP(i + n - 1, store=True)
op = self.ket.get_op(ops, i)
op = op.replace_labels(op_ax_p + op_ax_pstar, ax_p + ax_pstar)
C = self.ket.get_theta(i, n)
C = npc.tensordot(op, C, axes=[ax_pstar, ax_p]) # same labels
C = npc.tensordot(LP, C, axes=['vR', 'vL']) # axes_p + (vR*, vR)
C = npc.tensordot(C, RP, axes=['vR', 'vL']) # axes_p + (vR*, vL*)
C.ireplace_labels(['vR*', 'vL*'], ['vL', 'vR']) # back to original theta labels
theta_bra = self.bra.get_theta(i, n)
E.append(npc.inner(theta_bra, C, axes='labels', do_conj=True))
return np.real_if_close(np.array(E)) * self.bra.norm * self.ket.norm
def _contract_LP(self, i, LP):
"""Contract LP with the tensors on site `i` to form ``self._LP[i+1]``"""
LP = npc.tensordot(LP, self.ket.get_B(i, form='A'), axes=('vR', 'vL'))
axes = (self.ket._get_p_label('*') + ['vL*'], self.ket._p_label + ['vR*'])
# for a ususal MPS, axes = (['p*', 'vL*'], ['p', 'vR*'])
LP = npc.tensordot(self.bra.get_B(i, form='A').conj(), LP, axes=axes)
return LP # labels 'vR*', 'vR'
def _contract_RP(self, i, RP):
"""Contract RP with the tensors on site `i` to form ``self._RP[i-1]``"""
RP = npc.tensordot(self.ket.get_B(i, form='B'), RP, axes=('vR', 'vL'))
axes = (self.ket._get_p_label('*') + ['vR*'], self.ket._p_label + ['vL*'])
# for a ususal MPS, axes = (['p*', 'vR*'], ['p', 'vL*'])
RP = npc.tensordot(self.bra.get_B(i, form='B').conj(), RP, axes=axes)
return RP # labels 'vL', 'vL*'
def _to_valid_index(self, i):
"""Make sure `i` is a valid index (depending on `finite`)."""
if not self._finite:
return i % self.L
if i < 0:
i += self.L
if i >= self.L or i < 0:
raise ValueError("i = {0:d} out of bounds for MPSEnvironment".format(i))
return i
[docs]class TransferMatrix(sparse.NpcLinearOperator):
r"""Transfer matrix of two MPS (bra & ket).
For an iMPS in the thermodynamic limit, we often need to find the 'dominant `RP`' (and `LP`).
This mean nothing else than to take the transfer matrix of the unit cell and find the
(right/left) eigenvector with the largest (magnitude) eigenvalue, since it will dominate
:math:`(TM)^n RP` (or :math:`LP (TM)^n`) in the limit :math:`n \rightarrow \infty` - whatever
the initial `RP` is. This class provides exactly that functionality with :meth:`eigenvectors`.
Given two MPS, we define the transfer matrix as::
| ---M[i]---M[i+1]- ... --M[i+L]---
| | | |
| ---N[j]*--N[j+1]* ... --N[j+L]*--
Here the `M` denotes the matrices of the bra and `N` the ones of the ket, respectively.
To view it as a `matrix`, we combine the left and right indices to pipes::
| (vL.vL*) ->-TM->- (vR.vR*) acting on (vL.vL*) ->-RP
Note that we keep all M and N as copies.
Parameters
----------
bra : MPS
The MPS which is to be (complex) conjugated.
ket : MPS
The MPS which is not (complex) conjugated.
shift_bra : int
We start the `N` of the bra at site `shift_bra` (i.e. the `j` in the above network).
shift_ket : int | None
We start the `M` of the ket at site `shift_ket` (i.e. the `i` in the above network).
``None`` defaults to `shift_bra`.
transpose : bool
Wheter `self.matvec` acts on `RP` (``False``) or `LP` (``True``).
charge_sector : None | charges | ``0``
Selects the charge sector of the vector onto which the Linear operator acts.
``None`` stands for *all* sectors, ``0`` stands for the zero-charge sector.
Defaults to ``0``, i.e., *assumes* the dominant eigenvector is in charge sector 0.
form : ``'B' | 'A' | 'C' | 'G' | 'Th' | None`` | tuple(float, float)
In which canonical form we take the `M` and `N` matrices.
Attributes
----------
L : int
Number of physical sites involved in the transfer matrix, i.e. the least common multiple
of `bra.L` and `ket.L`.
shift_bra : int
We start the `N` of the bra at site `shift_bra`.
shift_ket : int | None
We start the `M` of the ket at site `shift_ket`. ``None`` defaults to `shift_bra`.
transpose : bool
Wheter `self.matvec` acts on `RP` (``True``) or `LP` (``False``).
qtotal : charges
Total charge of the transfer matrix (which is gauged away in matvec).
form : tuple(float, float) | None
In which canonical form (all of) the `M` and `N` matrices are.
flat_linop : :class:`~tenpy.linalg.sparse.FlatLinearOperator`
Class lifting :meth:`matvec` to ndarrays in order to use :func:`~tenpy.tools.math.speigs`.
pipe : :class:`~tenpy.linalg.charges.LegPipe`
Pipe corresponding to ``'(vL.vL*)'`` for ``transpose=False``
or to ``'(vR.vR*)'`` for ``transpose=True``.
label_split :
``['vL', 'vL*']`` if ``tranpose=False`` or ``['vR', 'vR*']`` if ``transpose=True``.
_bra_N : list of npc.Array
Complex conjugated matrices of the bra, transposed for fast `matvec`.
_ket_M : list of npc.Array
The matrices of the ket, transposed for fast `matvec`.
_contract_legs : int
Number of physical legs per site + 1.
"""
def __init__(self,
bra,
ket,
shift_bra=0,
shift_ket=None,
transpose=False,
charge_sector=0,
form='B'):
L = self.L = lcm(bra.L, ket.L)
if shift_ket is None:
shift_ket = shift_bra
self.shift_bra = shift_bra
self.shift_ket = shift_ket
self.transpose = transpose
if ket.chinfo != bra.chinfo:
raise ValueError("incompatible charges")
form = ket._to_valid_form(form)
p = ket._p_label # for ususal MPS just ['p']
assert p == bra._p_label
pstar = ket._get_p_label('*') # ['p*']
if not transpose: # right to left
label = '(vL.vL*)' # what we act on
label_split = ['vL', 'vL*']
M = self._ket_M = [
ket.get_B(i, form=form).itranspose(['vL'] + p + ['vR'])
for i in reversed(range(shift_ket, shift_ket + L))
]
N = self._bra_N = [
bra.get_B(i, form=form).conj().itranspose(pstar + ['vR*', 'vL*'])
for i in reversed(range(shift_bra, shift_bra + L))
]
pipe = npc.LegPipe([M[0].get_leg('vR'), N[0].get_leg('vR*')], qconj=-1).conj()
else: # left to right
label = '(vR*.vR)' # mathematically more natural
label_split = ['vR*', 'vR']
M = self._ket_M = [
ket.get_B(i, form=form).itranspose(['vL'] + p + ['vR'])
for i in range(shift_ket, shift_ket + L)
]
N = self._bra_N = [
bra.get_B(i, form=form).conj().itranspose(['vR*', 'vL*'] + pstar)
for i in range(shift_bra, shift_bra + L)
]
pipe = npc.LegPipe([N[0].get_leg('vL*'), M[0].get_leg('vL')], qconj=+1).conj()
dtype = np.promote_types(bra.dtype, ket.dtype)
self.pipe = pipe
self.label_split = label_split
self.flat_linop = sparse.FlatLinearOperator(self.matvec, pipe, dtype, charge_sector, label)
self.qtotal = bra.chinfo.make_valid(np.sum([B.qtotal for B in M + N], axis=0))
self._contract_legs = len(ket._p_label) + 1 # for a ususal MPS: 2
[docs] def matvec(self, vec):
"""Given `vec` as an npc.Array, apply the transfer matrix.
Parameters
----------
vec : :class:`~tenpy.linalg.np_conserved.Array`
Vector to act on with the transfermatrix.
If not `transposed`, `vec` is the right part `RP` of an environment,
with legs ``'(vL.vL*)'`` in a pipe or splitted.
If `transposed`, the left part `LP` of an environment with legs ``'(vR*.vR)'``.
Returns
-------
mat_vec : :class:`~tenpy.linalg.np_conserved.Array`
The tranfer matrix acted on `vec`, in the same form as given.
"""
pipe = None
if vec.rank == 1:
vec = vec.split_legs(0)
pipe = self.pipe
vec.itranspose(self.label_split) # ['vL', 'vL*'] or ['vR*', 'vR']
qtotal = vec.qtotal
legs = vec.legs
contract = self._contract_legs # number of physical legs per site + 1
# the actual work
if not self.transpose: # right to left
for N, M in zip(self._bra_N, self._ket_M):
vec = npc.tensordot(M, vec, axes=1) # axes=['vR', 'vL']
vec = npc.tensordot(vec, N, axes=contract) # [['p', 'vL*'], ['p*', 'vR*']]
else: # left to right
for N, M in zip(self._bra_N, self._ket_M):
vec = npc.tensordot(vec, M, axes=1) # axes=['vR', 'vL']
vec = npc.tensordot(N, vec, axes=contract) # [['vL*', 'p*'], ['vR*', 'p']])
if np.any(self.qtotal != 0):
# Hack: replace leg charges and qtotal -> effectively gauge `self.qtotal` away.
vec.qtotal = qtotal
vec.legs = legs
vec.test_sanity() # Should be fine, but who knows...
if pipe is not None:
vec = vec.combine_legs([0, 1], pipes=pipe)
return vec
[docs] def initial_guess(self, diag=1.):
"""Return a diagonal matrix as initial guess for the eigenvector.
Parameters
----------
diag : float | 1D ndarray
Should be ``1.`` for the identity or some singular values squared.
Returns
-------
mat : :class:`~tenpy.linalg.np_conserved.Array`
A 2D array with `diag` on the diagonal such that :meth:`matvec` can act on it.
"""
return npc.diag(diag, self.pipe.legs[0]).iset_leg_labels(self.label_split)
[docs] def eigenvectors(self,
num_ev=1,
max_num_ev=None,
max_tol=1.e-12,
which='LM',
v0=None,
**kwargs):
"""Find (dominant) eigenvector(s) of self using :mod:`scipy.sparse`.
If no charge_sector was selected, we look in *all* charge sectors.
Parameters
----------
num_ev : int
Number of eigenvalues/vectors to look for.
max_num_ev : int
:func:`scipy.sparse.linalg.speigs` somtimes raises a NoConvergenceError for small
`num_ev`, which might be avoided by increasing `num_ev`. As a work-around,
we try it again in the case of an error, just with larger `num_ev` up to `max_num_ev`.
``None`` defaults to ``num_ev + 2``.
max_tol : float
After the first `NoConvergenceError` we increase the `tol` argument to that value.
which : str
Which eigenvalues to look for, see `scipy.sparse.linalg.speigs`.
**kwargs :
Further keyword arguments are given to :func:`~tenpy.tools.math.speigs`.
Returns
-------
eta : 1D ndarray
The eigenvalues, sorted according to `which`.
w : list of :class:`~tenpy.linalg.np_conserved.Array`
The eigenvectors corresponding to `eta`, as npc.Array with LegPipe.
"""
if max_num_ev is None:
max_num_ev = num_ev + 2
flat_linop = self.flat_linop
if flat_linop.charge_sector is None:
# Try for all charge sectors
eta = []
A = []
for chsect in flat_linop.possible_charge_sectors:
flat_linop.charge_sector = chsect
eta_cs, A_cs = self.eigenvectors(num_ev, max_num_ev, max_tol, which, **kwargs)
eta.extend(eta_cs)
A.extend(A_cs)
flat_linop.charge_sector = None
else:
if v0 is not None:
kwargs['v0'] = self.flat_linop.npc_to_flat(v0)
# for given charge sector
for k in range(num_ev, max_num_ev + 1):
if k > num_ev:
warnings.warn("TransferMatrix: increased `num_ev` to " + str(k + 1))
try:
eta, A = speigs(flat_linop, k=k, which='LM', **kwargs)
A = np.real_if_close(A)
A = [flat_linop.flat_to_npc(A[:, j]) for j in range(A.shape[1])]
break
except scipy.sparse.linalg.eigen.arpack.ArpackNoConvergence:
if k == max_num_ev:
raise
# just retry with larger k and 'tol'
kwargs['tol'] = max(max_tol, kwargs.get('tol', 0))
# sort
perm = argsort(eta, which)
return np.array(eta)[perm], [A[j] for j in perm]
[docs]def build_initial_state(size, states, filling, mode='random', seed=None):
"""Build an "initial state" list.
Uses two iterables ('states' and 'filling') to determine how to fill the
state. The two lists should have the same length as every element in 'filling' gives the filling
fraction for the corresponding state in 'states'.
Example:
size = 6, states = [0, 1, 2], filling = [1./3, 2./3, 0.]
n_states = size * filling = [2, 4, 0]
==> Two sites will get state 0, 4 sites will get state 1, 0 sites will
get state 2.
.. todo ::
Make more general: it should be possible to specify states as strings.
Parameters
----------
size : int
length of state
states : iterable
Containing the possible local states
filling : iterable
Fraction of the total number of sites to get a certain state. If
infinite fractions (e.g. 1/3) are needed, one should supply a
fraction (1./3.)
mode : str | None
State filling pattern. Only 'random' is implemented
seed : int | None
Seed for random number generators
Returns
-------
initial_state (list) : the initial state
Raises
------
ValueError
If fractonal fillings are incommensurate with system size.
AssertionError
If the total filling is not equal to 1, or the length of `filling`
does not equal the length of `states`.
"""
random.seed(seed)
# Do some safety checks
assert sum(filling) == 1
assert len(states) == len(filling)
# Get number of sites for each local state
n_states = np.array(filling) * size
for num in n_states:
if ((num - round(num)) < 1e-12):
num = int(round(num))
if type(num) != int and not num.is_integer():
raise ValueError("Cannot create model of length {} with filling {}".format(
size, filling))
# Randomly assign local states
initial_state = [0] * size
all_sites = list(range(size)) # To avoid having two types on same site.
for state, fill in zip(states, filling):
sites = random.sample(set(all_sites),
int(fill * size)) # pick fill*size sites to put state
for site in sites:
initial_state[site] = state
all_sites.remove(site)
return initial_state