"""
Matrix Product Operators
************************
The matrix product operators organize the class :py:class:`mpo.MPO`. The
natural operator structure for MPS algorithms is a matrix product operator
(MPO). The exact diagonalization modules uses the same class to construct
Hamiltonians, effective Hamiltonians, and operators in Liouville space.
Operators can be constructed as MPOs using a set of rules which specify
how operators acting on local Hilbert spaces are combined to create an operator
acting on the entire many-body system. Using a canonical form for MPOs, several
operators constructed from these rules can be put together as a single
MPO :cite:`Wall2012`. OSMPS uses this essential idea of creating
operators from a set of predetermined rules to construct the Hamiltonian and
other operators appearing in an MPS simulation.
An operator acting on the many-body state, i.e. the MPO, is an object of the
:py:class:`mpo.MPO` class in MPSPyLib. For an instantiation of the MPO class,
named ``H`` in the following, is obtained as
.. code-block:: python
H = mps.MPO(Operators)
Once we have an MPO object, we add terms to it using the
:py:func:`Ops.AddMPOTerm` method. In particular, the rules are
defined in the module :py:mod:`mpoterm`.
.. automodule:: mpoterm
**Example**
As an example, a template MPO for the Hamiltonian of the hard-core boson model
with :math:`1/r^3` interactions,
.. math::
\hat{H} &= -t \sum_{\langle i,j \\rangle} \left(\hat{b}_i^{\\dagger} \hat{b}_j
+ \mathrm{h.c.} \\right)
+ U \sum_{i<j} \\frac{\hat{n}_i \hat{n}_j}{\left|j-i \\right|^3}
- \mu \sum_{i} \hat{n}_i \, ,
is provided by
.. code-block:: python
import MPSPyLib as mps
# Build operators
Operators = mps.BuildBoseOperators(1)
# Define Hamiltonian MPO
H = mps.MPO(Operators)
H.AddMPOTerm('site', 'nbtotal', hparam='mu', weight=-1.0)
H.AddMPOTerm('bond', ['bdagger', 'b'], hparam='t', weight=-1.0)
invrcube = lambda x: 1.0/(x*x*x)
H.AddMPOTerm('InfiniteFunction', ['nbtotal', 'nbtotal'], hparam='U',
weight=1.0, func=invrcube, L=1000, tol=1e-9)
We stress that no assumptions are made about the ``hparam`` parameters ``'t'``,
``'mu'``, or ``'U'`` when the MPO is built. The ``hparams`` are simply tags
parameterizing the operators which appear in the Hamiltonian. The actual
numerical values of :math:`t`, :math:`U`, and :math:`\mu` are specified as
part of the dictionary ``parameters`` discussed in Sec. :ref:`sec_parameters`.
These parameters can also be position-dependent, as demonstrated in the example
files in :ref:`sec_examples`.
One can print operators using the :py:func:`mpo.MPO.printMPO` method. As an
example, the interacting spinless fermion model
.. math::
\hat{H} &= - t \sum_{\langle i,j \\rangle} \left(\hat{f}_i^{\\dagger} \hat{f}_j
+ \mathrm{h.c.} \\right)
+ U\sum_{\langle i,j \\rangle }{\hat{n}_i \hat{n}_j} \, ,
is specified and printed as
.. code-block:: python
import MPSPyLib as mps
# Build operators
Operators = mps.BuildFermiOperators()
# Define Hamiltonian MPO
H = mps.MPO(Operators)
H.AddMPOTerm('bond', ['fdagger', 'f'], hparam='t', weight=-1.0, Phase=True)
H.AddMPOTerm('bond', ['nftotal', 'nftotal'], hparam='U', weight=1.0)
H.printMPO()
The result of this code is
:command:`-1.0 \* t \* \\sum_i fdaggerP_i f_{i+1} + 1.0 \* t \* \\sum_i fP_i`
:command:`fdaggerP_{i+1} + 1.0 \* U \* \\sum_i nftotal_i nftotal_{i+1}`
Note the use of the *phased* operators ``fdaggerP`` and ``fP`` according to
``Phase=True``.
As the ``'InfiniteFunction'`` rule constructs MPOs by fitting the decaying
function to a series of weighted exponentials, the result of ``printMPO``
will contain the data for the exponentials and not information about the
infinite function proper.
Finally, as complex MPOs with ``InfiniteFunction`` interactions can be
expensive to generate, one can save the data in an ``MPO`` object instance
using the :py:func:`mpo.WriteMPO` function, and then read it in later using
the :py:func:`mpo.ReadMPO` function. The syntax for these is
.. code-block:: python
mps.WriteMPO(H,filename)
H = mps.ReadMPO(filename)
where ``H`` is an object of the :py:class:`mpo.MPO` class and ``filename`` is
the name of the file where the MPO is written to/read from.
"""
import numpy as np
#import scipy.sparse as sp
from copy import deepcopy
import hashlib
#from .mpoterm import SiteTerm, BondTerm, ExpTerm, FiniteFunc, InfiniteFunc
#from .mpoterm import ProductTerm, MBStringTerm, TTerm
#from .mpoterm import SiteLindXX, MBSLind, MBSLindXY
#from .mpoterm import LindLSpec, LindSpec, LindInfFunc
from sys import version_info
if(version_info[0] < 3):
import cPickle
else:
import pickle as cPickle
from scipy import __version__ as _spv
_spv = list(map(int, _spv.split('.')[:2]))
interaction_picture = False
[docs]def set_interaction_picture(status):
global interaction_picture
interaction_picture = status
[docs]def WriteMPO(MPO, filename):
"""
Write MPO object instance to file filename.
**Arguments**
MPO : instance of :py:class:`mpo.MPO`
This is the MPO to be written.
filename : str
File where MPO information is to be written.
"""
output = open(filename, "wb")
cPickle.dump(MPO, output, cPickle.HIGHEST_PROTOCOL)
output.close()
[docs]def ReadMPO(filename):
"""
Read MPO object instance from file filename.
**Arguments**
filename : str
Filename where MPO is stored.
"""
output = open(filename, "rb")
MPO = cPickle.load(output)
output.close()
return MPO
[docs]class MPO:
"""
Initialize an MPO object.
**Arguments**
Operators : dict or instance of :py:class:`Ops.OperatorList`.
Contains the local operators used for measurements, symmetries,
and MPOs.
pbc : boolean
Contains the setting for open boundary conditions (``False``)
versus periodic boundary conditions (``True``). There are multiple
restrictions which methods can be used with PBC. In TN methods,
the only non-local terms supporting PBC are the bond rules and
the product rule. In ED, only bond rules access the pbc flag,
others might execute with OPC.
Default to False, i.e., open boundary conditions.
PostProcess : boolean, optional
If PostProcessing (``True``), infinite functions are not fitted
with exponentials. If ``False``, those are fitted.
**Variables**
TI : bool
Flag if MPO is translational invariant.
IntHparam : dictionary
Contains the mapping string identifier to integer for writting
to fortran interface.
terms : list
List of available MPO rules.
data : dictionary
Contains the class for each rule set with the string identifier
as key.
"""
def __init__(self, Operators, pbc=False, PostProcess=False):
self.Operators = Operators
self.PostProcess = PostProcess
self.TI = True
self.pbc = pbc
self.IntHparam = {}
self.terms = ['site', 'bond', 'exponential', 'product',
'FiniteFunction', 'MBString', 'tterm',
'lind1', 'MBSLind', 'MBSLXY', 'lindexp',
'InfiniteFunction', 'LindInfFunc']
self.data = {}
self.data['site'] = SiteTerm(Operators)
self.data['bond'] = BondTerm(Operators)
self.data['product'] = ProductTerm(Operators)
self.data['exponential'] = ExpTerm(Operators)
self.data['FiniteFunction'] = FiniteFunc(Operators)
self.data['MBString'] = MBStringTerm(Operators)
self.data['lind1'] = SiteLindXX(Operators)
self.data['tterm'] = TTerm(Operators)
self.data['MBSLind'] = MBSLind(Operators)
self.data['MBSLXY'] = MBSLindXY(Operators)
self.data['lindexp'] = ExpTerm(Operators)
self.data['InfiniteFunction'] = InfiniteFunc(Operators,
self.data['exponential'],
PostProcess)
self.data['LindInfFunc'] = LindInfFunc(Operators, self.data['lindexp'],
PostProcess)
self.data['lindlspec'] = LindLSpec(Operators)
self.data['lindspec'] = LindSpec(Operators)
[docs] def has_symmetry(self, Generator, has_cmplx, isz2, eps=1e-8):
"""
Check if the specific generator commutes with the terms in the
Hamiltonian. Returns ``True`` (``False``) if Hamiltonian respects
(does not fulfill) symmetry of the generator.
**Arguments**
Generator : matrix (np.array)
Generators of the symmetry should commute with the Hamiltonian.
Do not pass the generators as list, but loop over them and call
this function for each of them.
has_cmplx : Bool
boolean if any operator contains complex values.
isz2 : bool
Check for Z2 symmetry (True), otherwise (False) check for U1.
eps : float, optional
Accept non-zero commutators up to `eps`, default to
:math:`10^{-8}`.
**Details**
The checks for the symmetry follow the following scheme.
* For the site terms every operator :math:`O` is checked that it
fulfills :math:`[G, O] = 0`, where :math:`G` is the generator.
* For the bond terms the commutator :math:`[G2, H_{i,j}] = 0` is
considered with :math:`G2 = G \otimes 1 + 1 \otimes G` and
:math:`H_{i,j}` is the sum over all bond terms and includes therefore
the hermitian conjugate terms of every contribution as well.
* The product terms are not tested. The corresponding matrix
on the complete Hilbert space necessary for the operator could
exceed computational resources.
* Finite functions consist of two term and are tested according to
the procedure for bond terms. Only the nearest neighbor terms
are tested
* The exponential rule consists of two term and are tested according
to the procedure for bond terms. Only the nearest neighbor terms
are tested.
* Infinite functions consist of two term and are fitted to exponentials
anyway. So see for exponentials.
* MBString terms might be to big to be checked on the corresponding
Hilbert space.
The Fermi phase is accounted for in the MPO terms except for the
intermediate sites. For the generator the Fermi phase operators are
not necessary assuming that the generator consists of
:math:`b^{\dagger} b` with canceling terms :math:`\sigma^{z}` during
the Jordan-Wigner transformation.
"""
# Define function to calculate the commutators
comm_sum = lambda x, y : np.sum(np.abs(np.dot(x, y) - np.dot(y, x)))
# Check site terms
# ----------------
try:
for ii in range(len(self.data['site']['Op'])):
if(comm_sum(Generator,
self.Operators[self.data['site']['Op'][ii]])
> eps):
return False
except KeyError:
pass
# Check two-site terms (only nearest neighbors)
# ---------------------------------------------
#
# * bond terms
# * finite functions
# * exponentials
keys = ['bond', 'FiniteFunction', 'exponential']
# Build generator on two sites for all the following cases
Gen2 = np.kron(Generator, np.eye(Generator.shape[0])) + \
np.kron(np.eye(Generator.shape[0]), Generator)
if(isz2): Gen2 = np.mod(Gen2, 2)
for key in keys:
if(len(self.data[key]) > 0):
# Build sum of all operators (to include h.c.)
dim = self.Operators[self.data[key]['Op'][0][0]].shape
if(has_cmplx):
Ops = np.zeros((dim[0]**2, dim[1]**2), dtype=np.complex128)
else:
Ops = np.zeros((dim[0]**2, dim[1]**2))
for ii in range(len(self.data[key]['Op'])):
name1 = self.data[key]['Op'][ii][0]
name2 = self.data[key]['Op'][ii][1]
weight = self.data[key]['weight'][ii]
Ops += weight * np.kron(self.Operators[name1],
self.Operators[name2])
if(comm_sum(Gen2, Ops) > eps):
return False
# Check product terms
# -------------------
#
# skipped (space might get to big)
# Check MBString terms
# --------------------
#
# skipped (space might get to big)
return True
[docs] def get_param_set(self):
"""
Get a set with all the keys appearing in the definition of the MPO
for statics (hashing).
"""
params = set()
for key in self.terms:
if(key in ['lind1']): continue
params |= set(self.data[key]['hp'])
return list(params)
[docs] def get_hash(self):
"""
Get a hash function identifying the whole object.
"""
hstr = str(self.TI)
for key in self.terms:
hstr += self.data[key].get_hash()
hstr += self.data['InfiniteFunction'].get_hash()
hh = hashlib.sha512()
hh.update(str.encode(hstr))
return hh.hexdigest()
[docs] def printMPO(self):
"""
Print out the MPO.
"""
tstring = ''
for key in self.terms:
tstring += self.data[key].print_term()
print(tstring)
[docs] def AddMPOTerm(self, termtype, terms, **kwargs):
"""
Add a term to an MPO object
**Arguments**
termtype : str
Identifies which MPO term should be added. The following options
are available: `site`, `bond`, `FiniteFunction`,
`exponential`, `product`, `MBString`, `InfiniteFunction`, and
``tterm``.
terms : str, list of strings
This arguments varies from MPOTerm to MPO term. You can find the
explicit description for each term in the corresponding
description
* :py:class:`mpoterm.SiteTerm`
* :py:class:`mpoterm.BondTerm`
* :py:class:`mpoterm.ExpTerm`
* :py:class:`mpoterm.FiniteFunc`
* :py:class:`mpoterm.InfiniteFunc`
* :py:class:`mpoterm.ProductTerm`
* :py:class:`mpoterm.MBStringTerm`
* :py:class:`mpoterm.TTerm`
* :py:class:`mpoterm.SiteLindXX` (exact diagonalization)
``**kwargs`` : depends
The optional keyword argument depend as well on the type of the
MPO term. See links to descriptions in previous arguments.
"""
try:
self.data[termtype].AddTerm(terms, **kwargs)
except KeyError:
raise Exception("Unknown MPO term type " + str(termType)
+ " passed to AddMPOTerm!")
# -------------------------- Other functions -------------------------------
[docs] def getHparams(self):
"""
Get all the unique Hamiltonian parameters which appear in the MPO
definition.
"""
# NOTE: InfiniteFunctions are actually exponentials for MPS
# empty set
hparams = set()
# add all unique elements of any type to this set
for key in self.terms:
hparams |= set(self.data[key]['hp'])
return hparams
[docs] def IndexHParams(self):
"""
Index the hamiltonian parameters using integers instead of
human-readable keys.
"""
# Get list of hparams and put 'one' at the end
hparams = list(self.getHparams())
last = len(hparams)
if('one' in hparams):
last -= 1
idx = hparams.index('one')
if(idx != len(hparams) - 1):
hparams[idx] = hparams[-1]
hparams[-1] = 'one'
# Give Hamiltonian Parameters unique integer keys
self.IntHparam = {}
count = 0
for param in hparams[:last]:
self.IntHparam[count] = param
self.IntHparam[param] = count + 1
count += 1
self.IntHparam['num'] = count
return
[docs] def write(self, Openfile, Parameters, IntOperators, IntHparam):
"""
Write the MPO rules to a Fortran-readable file.
**Arguments**
Openfile : filehandle
open file where the parameters should be written.
Parameters : dictionary
dictionary with parameters for one simulation.
IntOperators : dictionary
mapping the string identifiers for the operators to integer
identifiers used within fortran.
IntHparam : dictionary
mapping the string identifiers for Hamiltonian parameters to
integer identifiers used within fortran.
"""
# Check if the MPO is TI
# ----------------------
TI = True
for key in self.terms[:11]:
TI = TI and self.data[key].check_ti(Parameters)
if((not TI) and (Parameters['simtype'] == 'Infinite')):
raise Exception("Site-dependent Hamiltonians not currently "
+ "supported for infinite-size systems!")
TI = 'T\n' if (TI) else 'F\n'
Openfile.write(TI)
pbc_str = 'T\n' if (self.pbc) else 'F\n'
Openfile.write(pbc_str)
# Write out terms for each rule
for key in self.terms[:11]:
self.data[key].write(Openfile, IntHparam, IntOperators, Parameters)
return
# -------------------- Hamiltonian on complete space -----------------------
[docs] def build_hamiltonian(self, param, SymmSec, sparse, maxdim,
quenchtime=None):
"""
Build the Hamiltonian for the simulation.
**Arguments**
param : dictionary with simulation parameters
settings of the whole simulation.
SymmSec : instance of SymmmetrySector or None
No meaning here, needed to keep calls equal to sparse equivalent
sparse : bool
flag if sparse matrix should be used for Hamiltonian.
maxdim : int
limits the dimension of the complete Hilbert space.
quenchtime : list / tuple, optional
If present first entries specifies the number of the quench and
second entry the actual time in the whole time evolution. Only
necessary for Hamiltonian of dynamics.
Default to None.
"""
dtype = np.complex128 if(self.Operators.has_complex) else np.float64
if(not SymmSec is None):
# Case with symmetries
return self._build_hamiltonian_symm(param, SymmSec, sparse, maxdim,
quenchtime)
return self._build_hamiltonian_nosymm(param, sparse, maxdim,
quenchtime)
[docs] def hparam(self, hpstr, param, quenchtime):
"""
Return the actual parameter/coupling for a term in the
Hamiltonian.
**Arguments**
hpstr : str
name of the operator
param : dict
setup for the simulation. Contains information about couplings
and dynamics.
quenchtime : None or tuple of int and float
For statics pass None. For dynamics pass integer ii for the
ii-th quench and the actual time t as float in a tuple as
(ii, t).
"""
if(hpstr == 'one'):
hp = 1.0
elif(quenchtime is None):
hp = param[hpstr]
else:
quench = quenchtime[0]
time = quenchtime[1]
if(hpstr in param['Quenches'].data[quench]['Hparams']):
# Time-dependent hparam
idx = param['Quenches'].data[quench]['Hparams'].index(hpstr)
hp = param['Quenches'].data[quench]['funcs'][idx](time)
else:
# Time-independent hparam (use initial value
hp = param[hpstr]
if(not hasattr(hp, '__len__')):
hp = np.ones((param['L'])) * hp
else:
hp = np.array(deepcopy(hp))
return hp
[docs] def _col_ind_dic(self, Op):
"""
Create a dictionary taking a row index as key and returning the
corresponding non-zero elements of the operators as tuple
(column index, value).
**Arguments**
Op : numpy 2d array
Operators to be transfered into this format.
"""
dic = {}
for i1 in range(Op.shape[0]):
for i2 in range(Op.shape[1]):
if(Op[i1, i2] != 0.0):
try:
dic[i1].append((i2, Op[i1, i2]))
except KeyError:
dic[i1] = [(i2, Op[i1, i2])]
return dic
[docs] def _build_hamiltonian_nosymm(self, param, sparse, maxdim,
quenchtime):
"""
Build the Hamiltonian on a complete Hilbert space as sparse/dense
matrix. For arguments description look into
:py:func:`mpo.MPO.build_hamiltonian`.
"""
# Get shortcuts to system size and local dimension of Hilbert space
ll = param['L']
ld = self.Operators['I'].shape[0]
dim = ld**ll
# Data type (always complex if not purely dissipative)
dtype = np.complex128 if(self.Operators.has_complex) else np.float64
# Initialize matrix in dependence of sparse/dense
if(sparse):
spm = sp.csc_matrix if(not quenchtime is None) else sp.csr_matrix
Ham = spm((dim, dim), dtype=dtype)
else:
Ham = np.zeros((dim, dim), dtype=dtype)
if((not maxdim is None) and (dim > maxdim)):
raise Exception('Dimension exceeds maximal dimension ' +
'of the Hilbert space specified.')
keylist = ['site', 'bond', 'exponential', 'MBString', 'tterm',
'InfiniteFunction']
for key in keylist:
if(len(self.data[key]) == 0): continue
for elem in self.data[key].build_nosymm(param, sparse, quenchtime,
ll, ld, self.pbc):
Ham += elem
keylist = ['product', 'FiniteFunction']
for key in keylist:
if(len(self.data[key]) > 0):
raise NotImplementedError('Term ' + key +
+ ' is not implemented for ED!')
return Ham
[docs] def _build_hamiltonian_symm(self, param, SymmSec, sparse, maxdim,
quenchtime, lil=False):
"""
Build a sparse or dense Hamiltonian directly in the symmetry-adapted
subspace. For argument description see
:py:func:`mpo.MPO.build_hamiltonian`.
"""
ll = param['L']
spdim = len(SymmSec)
dtype = np.complex128 if(self.Operators.has_complex) else np.float64
if(sparse):
SpHam = sp.lil_matrix((spdim, spdim), dtype=dtype)
else:
SpHam = np.zeros((spdim, spdim), dtype=dtype)
keylist = ['site', 'bond', 'exponential', 'MBString', 'InfiniteFunction']
for key in keylist:
if(len(self.data[key]) == 0): continue
for k1, k2, elem in self.data[key].build_symm(param, SymmSec,
quenchtime, ll,
self.pbc):
SpHam[k1, k2] += elem
keylist = ['product', 'FiniteFunction', 'tterm']
for key in keylist:
if(len(self.data[key]) > 0):
raise NotImplementedError('Term ' + key +
+ ' is not implemented for ED!')
if(sparse and lil): return SpHam
if(sparse and (not quenchtime is None)): return sp.csc_matrix(SpHam)
if(sparse): sp.csr_matrix(SpHam)
return SpHam
# --------------------- Build effective Hamiltonian ------------------------
[docs] def build_effhamiltonian(self, param, SymmSec, sparse, maxdim, quenchtime):
"""
Build the effective Hamiltonian representing the hermitian part of the
Lindblad master equation.
**Arguments**
param : dictionary with simulation parameters
settings of the whole simulation.
SymmSec : instance of SymmmetrySector or None
No meaning here, needed to keep calls equal to sparse equivalent
sparse : bool
flag if sparse matrix should be used for Hamiltonian.
maxdim : int
limits the dimension of the complete Hilbert space.
quenchtime : list / tuple
If present first entries specifies the number of the quench and
second entry the actual time in the whole time evolution. The
effective Hamiltonian can only be used in the dynamics, so term
has to be present.
"""
if(not SymmSec is None):
# Case with symmetries
return self._build_effhamiltonian_symm(param, SymmSec, sparse,
maxdim, quenchtime)
return self._build_effhamiltonian_nosymm(param, sparse, maxdim,
quenchtime)
[docs] def _build_effhamiltonian_nosymm(self, param, sparse, maxdim,
quenchtime):
"""
Build the effective Hamiltonian describing the hermitian part of
the Lindblad master equation. For description of arguments see
:py:func:`mpo.MPO.build_effhamiltonian`.
"""
# Get shortcuts
ll = param['L']
ld = self.Operators['I'].shape[0]
dim = ld**ll
# Define functions in dependence of sparse/dense; initialize matrix
if(sparse):
fmt = 'csc' if(not quenchtime is None) else 'csr'
if(fmt == 'csc'):
tocmplx = lambda Mat : sp.csc_matrix(Mat, dtype=np.complex128)
else:
tocmplx = lambda Mat : sp.csr_matrix(Mat, dtype=np.complex128)
else:
tocmplx = lambda Mat : np.array(Mat, dtype=np.complex128)
# Get Hamiltonian - convert to complex array
Ham = self._build_hamiltonian_nosymm(param, sparse, maxdim, quenchtime)
if(Ham.dtype != np.complex128):
Ham = tocmplx(Ham)
if(self.data['lind1']):
# Loop over local Lindblads
for elem in self.data['lind1'].build_nosymm(param, sparse,
quenchtime, ll, ld,
self.pbc):
Ham += 1j * elem
return Ham
[docs] def _build_effhamiltonian_symm(self, param, SymmSec, sparse, maxdim,
quenchtime, lil=False):
"""
Build the effective Hamiltonian decribing the hermitian part of the
Lindblad master equation. For description of arguments see
:py:func:`mpo.MPO.build_effhamiltonian`.
"""
# Get shortcuts
ll = param['L']
ld = self.Operators['I'].shape[0]
spdim = len(SymmSec)
# Get Hamiltonian and convert to complex matrix
dtype = np.complex128
SpHam = self._build_hamiltonian_symm(param, SymmSec, sparse, maxdim,
quenchtime, lil=True)
if(sparse and (SpHam.dtype != dtype)):
SpHam = sp.lil_matrix(SpHam, dtype=dtype)
elif(SpHam.dtype != dtype):
SpHam = np.array(SpHam, dtype=dtype)
# Loop over effective Lindblad terms
if(len(self.data['lind1']) > 0):
for k1, k2, elem in self.data['lind1'].build_symm(param, SymmSec,
quenchtime, ll):
SpHam[k1, k2] += 1j * elem
if(sparse and lil): return SpHam
if(sparse): return sp.csc_matrix(SpHam)
return SpHam
# ---------------------- Build Liouville operators -------------------------
[docs] def build_liouville(self, param, SymmSec, sparse, maxdim, quenchtime=None):
"""
Build the Liouville operator for the simulation.
**Arguments**
param : dictionary with simulation parameters
settings of the whole simulation.
SymmSec : instance of SymmmetrySector or None
No meaning here, needed to keep calls equal to sparse equivalent
sparse : bool
flag if sparse matrix should be used for Hamiltonian.
maxdim : int
limits the dimension of the complete Hilbert space.
quenchtime : list / tuple, optional
If present first entries specifies the number of the quench and
second entry the actual time in the whole time evolution. Only
necessary for Hamiltonian of dynamics.
Default to None.
"""
if(not SymmSec is None):
# Case with symmetries
return self._build_liouville_symm(param, SymmSec, sparse, maxdim,
quenchtime)
return self._build_liouville_nosymm(param, sparse, maxdim, quenchtime)
[docs] def _build_liouville_nosymm(self, param, sparse, maxdim, quenchtime):
"""
Build the Liouville operator without any symmetry conservation. For
arguments look into :py:func:`mpo.MPO.build_liouville`.
"""
# Get shortcuts to system size and local dimension of Hilbert space
ll = param['L']
ld = self.Operators['I'].shape[0]
dim = ld**ll
# Data type (always complex if not purely dissipative)
dtype = np.complex128# if(self.Operators.has_complex) else np.float64
# Define functions in dependence of sparse/dense; initialize matrix
if(sparse):
fmt = 'csc' if(not quenchtime is None) else 'csr'
kron = lambda x, y: sp.kron(x, y, format=fmt)
if((_spv[0] < 1) and (_spv[1] < 12)):
eye = lambda dim : sp.eye(dim, dim)
else:
eye = sp.eye
else:
kron = np.kron
eye = np.eye
Ham = self._build_hamiltonian_nosymm(param, sparse, maxdim,
quenchtime)
Liou = kron(- 1j * Ham, eye(dim)) \
+ kron(eye(dim), 1j * Ham.transpose())
if(len(self.data['lindspec']) > 0):
# Get Lindblad operators from spectral decomposition
global interaction_picture
if(interaction_picture): Liou[:, :] = 0.0
for elem in self.data['lindspec'].lbuild_nosymm(param, sparse,
quenchtime, ll, ld,
Ham):
Liou += elem
return Liou
if(len(self.data['lindlspec']) > 0):
# Get Lindblad operators from spectral decomposition of a local
# block
for elem in self.data['lindlspec'].lbuild_nosymm(param, sparse,
quenchtime, ll, ld,
self):
Liou += elem
return Liou
if(len(self.data['lind1']) > 0):
for elem in self.data['lind1'].lbuild_nosymm(param, sparse,
quenchtime, ll, ld,
self.pbc):
Liou += elem
if(len(self.data['MBSLind']) > 0):
for elem in self.data['MBSLind'].lbuild_nosymm(param, sparse,
quenchtime, ll, ld,
self.pbc):
Liou += elem
if(len(self.data['MBSLXY']) > 0):
for elem in self.data['MBSLXY'].lbuild_nosymm(param, sparse,
quenchtime, ll, ld,
self.pbc):
Liou += elem
return Liou
[docs] def _build_liouville_symm(self, param, SymmSec, sparse, maxdim, quenchtime):
"""
Build propagator in Liouville space with symmetry. For
arguments look into :py:func:`mpo.MPO.build_liouville`.
"""
# Shortcuts
ll = param['L']
dim = len(SymmSec)
# Data type (always complex if not purely dissipative)
dtype = np.complex128
# Define functions in dependence of sparse/dense; initial matrix
if(sparse):
kron = lambda x, y : sp.kron(x, y, format='csr')
else:
kron = lambda x, y : np.kron(x, y)
eye = sp.eye if(sparse) else np.eye
Ham = self._build_hamiltonian_symm(param, SymmSec, sparse, maxdim,
quenchtime)
Liou = kron(Ham, -1j * eye(dim)) + kron(1j * eye(dim), Ham.transpose())
if(sparse): Liou = Liou.tolil()
if(len(self.data['lindspec']) > 0):
# Get Lindblad operators from spectral decomposition
global interaction_picture
if(interaction_picture): Liou[:, :] = 0.0
for elem in self.data['lindspec'].lbuild_symm(param, SymmSec,
quenchtime, ll, Ham):
Liou += elem
return Liou
if(self.data['lind1']):
# Loop over local terms
for k1, k2, elem in self.data['lind1'].lbuild_symm(param, SymmSec,
quenchtime, ll,
self.pbc):
Liou[k1, k2] += elem
if(self.data['MBSLind']):
for k1, k2, elem in self.data['MBSLind'].lbuild_symm(param, SymmSec,
quenchtime, ll,
self.pbc):
Liou[k1, k2] += elem
if(self.data['MBSLXY']):
for k1, k2, elem in self.data['MBSLXY'].lbuild_symm(param, SymmSec,
quenchtime, ll,
self.pbc):
Liou[k1, k2] += elem
if(sparse and (not quenchtime is None)): return sp.csc_matrix(Liou)
if(sparse): sp.csr_matrix(Liou)
return Liou