"""
MPO Terms to build rule sets
****************************
The following MPO terms are available at present
* ``site`` rule (:py:class:`mpoterm.SiteTerm`). Example:
.. code-block:: python
H.AddMPOTerm('site', 'sz', hparam='h_z', weight=1.0)
generates :math:`\hat{H}=h_z\sum_{i=1}^{L}\hat{S}^z_i`, where
:math:`\hat{S}^z =` ``Operators['sz']}`` and ``Operators`` was already
passed when creating the instance of the MPO. The second argument contains the
relevant keys of the operators to be used. Because only the names of the
operators are used, this routine handles ordinary and symmetry-adapted
operators on the same footing.
* ``bond`` rule (:py:class:`mpoterm.BondTerm`). Example:
.. code-block:: python
H.AddMPOTerm('bond', ['splus','sminus'], hparam='J_xy', weight=0.5)
generates :math:`\hat{H} = \\frac{J_{xy}}{2} \sum_{\langle i,j \\rangle}`
:math:`\left(\hat{S}^{+}_{i} \hat{S}^{-}_{j} + \mathrm{h.c.} \\right)`,
where :math:`\langle i,j\\rangle` denotes all nearest neighbor pairs :math:`i`
and :math:`j`,
.. code-block:: python
H.AddMPOTerm('bond', ['sz', 'sz'], hparam='J_z', weight=1.0)
generates :math:`\hat{H}=J_{z}\sum_{\langle i,j\\rangle}\hat{S}^z_i\hat{S}^z_j`,
and
.. code-block:: python
H.AddMPOTerm('bond', ['fdagger', 'f'], hparam='t', weight=-1.0, Phase=True)
generates :math:`\hat{H} = -t \sum_{\langle i,j \\rangle}`
:math:`\left(\hat{f}_i^{\dagger} \hat{f}_j + \mathrm{h.c.} \\right)`,
where :math:`\hat{f}` is a fermionic destruction operator. That is, terms
with Fermi phases are cared for with the optional ``Phase`` argument
and Hermiticity is automatically enforced by this routine. The Fermi phase
is added through a Jordan-Wigner string of the operators in ``Operator``
with key ``FermiPhase``. The functions :py:func:`Ops.BuildFermiOperators`
and :py:func:`Ops.BuildFermiBoseOperators` discussed in
Sec. :ref:`sec_deflocalops` both create the Fermi phase operator appropriate
for this use automatically. Notice that the order of the operators matters
due to the Jordan-Wigner transformation. Adding ``['f', 'fdagger']`` does not
result in the same fermionic Hamiltonian. A warning will be printed for obvious
cases.
* ``exponential`` rule set (:py:class:`mpoterm.ExpTerm`). Example
.. code-block:: python
H.AddMPOTerm('exponential', ['sz', 'sz'], hparam='J_z', decayparam=a, weight=1.0)
generates :math:`\hat{H} =J_z \sum_{i<j} a^{j-i+1} \hat{S}^z_i \hat{S}^z_j`.
This rule also enforces Hermiticity and cares for Fermi phases as in the
bond rule case.
* ``FiniteFunction`` rule (:py:class:`mpoterm.FiniteFunc`). Example:
.. code-block:: python
f = []
for ii in range(6):
f.append(1.0 / (ii + 1.0)**3)
H.AddMPOTerm('FiniteFunction', ['sz', 'sz'], f=f, hparam='J_z', weight=-1.0)
generates :math:`\hat{H} = -J_z \sum_{1 \le i,j \le 6}`
:math:`\\frac{1}{\left|j-i \\right|^3} \hat{S}^z_i \hat{S}^z_j`,
where the summation is over all :math:`i` and :math:`j` pairs separated by
at least :math:`1` and at most :math:`6` lattice spacings. The range of the
potential is given by the number of elements in the list ``f``. This rule
enforces Hermiticity and cares for Fermi phases as in the bond rule case.
* ``InfiniteFunction`` rule set (:py:class:`mpoterm.InfiniteFunc`). We point
out that this term is always fitted with a set of :py:class:`mpoterm.ExpTerm`
for MPS simulations. Examples:
.. code-block:: python
invrcube = lambda x: 1.0/(x*x*x)
H.AddMPOTerm('InfiniteFunction', ['sz', 'sz'], hparam='J_z', weight=1.0,
func=invrcube, L=1000, tol=1e-9)
generates
:math:`\hat{H}=\sum_{i<j}\\frac{1}{\left|j-i\\right|^3}\hat{S}^z_i\hat{S}^z_j`,
where the functional form is valid to a distance ``L`` with a residual of at
most ``tol``. Similarly,
.. code-block:: python
H.AddMPOTerm('InfiniteFunction', ['sz', 'sz'], hparam='J_z', weight=1.0,
x=x, y=f, tol=1e-9)
generates :math:`\hat{H}=\sum_{i<j}f\left(j-i\\right)\hat{S}^z_i\hat{S}^z_j`,
where the functional form :math:`f\left(x\\right)` is determined by the array
of values ``f`` and the array of evaluation points ``x``. The distance of
validity ``L`` is determined from the array ``x``. This rule applies only
to monotonically decreasing functions ``func`` or ``f``, and so outside the
range of validity ``L`` the corrections are also monotonically decreasing.
This rule enforces Hermiticity and cares for Fermi phases as in the bond rule
case.
* ``product`` rule (:py:class:`mpoterm.ProductTerm`). Example:
.. code-block:: python
H.AddMPOTerm('product', 'sz', hparam='h_p', weight=-1.0)
generates :math:`\hat{H}=-h_p\prod_{i=1}^{L} \hat{S}^z_i`. The operator used
in the product must be Hermitian.
* ``MBString`` rule set(:py:class:`mpoterm.MBStringTerm`). See function
description for usage.
* :py:class:`mpoterm.TTerm`. See function description for usage.
* :py:class:`mpoterm.SiteLindXX`. See function description for usage.
* :py:class:`mpoterm.LindExp`. See function description for usage.
* :py:class:`mpoterm.LindInfFunc`. See function description for usage.
* :py:class:`mpoterm.MBSLind`. See function description for usage.
(Can act as nearest neighbor Lindblad term.)
In all of the above expressions, ``H`` is an object of the ``MPO`` class.
The strings passed to the argument ``hparams`` represent tags for variable
Hamiltonian parameters, which allow for easy batching and parallelization
(Sec. :ref:`sec_parallel`), and also facilitate dynamical processes
The argument ``hparams`` is optional, and defaults to the constant
value 1. The optional argument ``weight`` specifies a constant factor that
scales the operator and defaults to the value 1, too.
**Developer's Guide for Rule Sets**
Module contains a class for each possible term in the MPO rule sets. Each term
is derived from a base class.If you plan to add additional MPO terms to OSMPS,
you have to modify the following files. In the case of a new Lindblad operator,
these are
* mpoterm.py: add class for the new term
* mpo.py: add new term to the MPO structure and enable writing.
* MPOOps_types_templates.f90: define derived type and add to MPORuleSet
* MPOOps_include.f90: add read/destroy methods for RuleSets
* MPOOps_include.f90: effective Hamiltonian for 2 sites
* MPOOps_include.f90: Liouville operator for 2 sites
* MPOOps_include.f90: effective Hamiltonian represented as MPO.
* MPOOps_include.f90: Liouville operator represented as MPO.
* TimeevolutionOps_include.f90: Add overlaps to QT functionality.
"""
import numpy as np
#import numpy.linalg as la
#import scipy.sparse as sp
#import scipy.linalg as sla
#from scipy.optimize import leastsq
from copy import deepcopy
import hashlib
from itertools import product as Cartesian
#from MPSPyLib.EDLib import CouplSpinInBEC, CouplSpinInEM3D
from scipy import __version__ as _spv
_spv = list(map(int, _spv.split('.')[:2]))
# Original herm check: add if (not A^dagger == A) and (not B^dagger == B)?
# Decomposition: whole operator has to be hermitian
# InfiniteFunc: x, y, func, L check could be better
# Where is TI used as class variable in MPO
[docs]def _DecomposeOperator(Operator):
"""
Performs a singular value decomposition on a bond operator
:math:`O_{[i j], [i' j']}` to find a Frobenius orthogonalized set
of operators :math:`O^L` and :math:`O^R` such that
:math:`O_{[i j], [i' j']} = \sum_a (O^L_a)_{i i'} (O^R_a)_{j j'}`.
This allows us to only consider bond terms that are tensor products.
**Arguments**
Operator : 2d numpy array
Operator is represented as np 2d matrix of dimension
:math:`d^2 \\times d^2` where :math:`d` is the local dimenison
of one site in the Hilbert space.
"""
d = int(np.sqrt(Operator.shape[0]))
reshapedOp = np.transpose(np.reshape(Operator, [d, d, d, d]), [0, 2, 1, 3])
reshapedOp = np.reshape(reshapedOp, [d**2, d**2])
U, S, Vh = sla.svd(reshapedOp)
eps = np.finfo(np.double).eps
dn1 = np.dot(S,S)
dn2 = 0.0
count = 0
for si in S:
count += 1
dn2 += si**2
if(1 - dn2 / dn1 < eps): break
OpList = []
for i in range(count):
OL = np.reshape(U[:, i] * np.sqrt(S[i]), [d, d])
OR = np.reshape(Vh[i, :] * np.sqrt(S[i]), [d, d])
OpList.append([OL, OR])
return OpList
[docs]def DecomposeOperator(Operator, gen):
"""
Performs a singular value decomposition on a bond operator
:math:`O_{[i j], [i' j']}` to find a Frobenius orthogonalized set
of operators :math:`O^L` and :math:`O^R` such that
:math:`O_{[i j], [i' j']} = \sum_a (O^L_a)_{i i'} (O^R_a)_{j j'}`.
This allows us to only consider bond terms that are tensor products.
**Arguments**
Operator : 2d numpy array
Operator is represented as np 2d matrix of dimension
:math:`d^2 \\times d^2` where :math:`d` is the local dimenison
of one site in the Hilbert space.
gen : list of np 1d-arrays or ``None``
Contains the diagonal entries of the diagonal symmetry generator
for a two-site operator for each symmetry. If there are no
symmetries, pass ``None``.
"""
if(gen is None):
# Case without symmetry generator
# -------------------------------
return _DecomposeOperator(Operator)
# Case with symmetry generator
# ----------------------------
qnumbers = []
for ii in range(gen[0].shape[0]):
genii = []
for elem in gen:
genii.append(elem[ii])
qnumbers.append(tuple(genii))
# That are the unique quantum numbers
qnumbers_unique = set(qnumbers)
# Collect the indices with the same quantum numbers
qidx = []
for elem in qnumbers_unique:
qidx.append([])
for ii in range(gen[0].shape[0]):
if(elem == qnumbers[ii]):
qidx[-1].append(ii)
# Build operator for each quantum number and collect in list
OpList = []
for ii in range(len(qidx)):
Oq = np.zeros(Operator.shape, dtype=Operator.dtype)
for jj in range(len(qidx[ii])):
Oq[qidx[ii][jj], qidx[ii]] = Operator[qidx[ii][jj], qidx[ii]]
Qlist = _DecomposeOperator(Oq)
for jj in range(len(Qlist)):
Qlist[jj][0][np.abs(Qlist[jj][0]) < 1e-15] = 0.0
Qlist[jj][1][np.abs(Qlist[jj][1]) < 1e-15] = 0.0
OpList.append([Qlist[jj][0], Qlist[jj][1]])
return OpList
[docs]def fittoexpSum(x, y, maxnterms, tol):
"""
Fit a sum of exponentials :math:`\sum_n a_n b_n^{x-1}` to the function
across the range [1:L].
**Arguments**
x : 1d array
represents the points x where the function f(x) is evaluated
and fitted.
y : 1d array
Actual values of f(x) which should be represented by the sum
of exponentials.
int : maxnterms
the maximum number of exponentials allowed
float : tol
the tolerance used to obtain the actual number of terms.
"""
# function evaluated at x-points
maxnterms = min(maxnterms, (len(x) // 2))
fail = True
for n in range(1, maxnterms + 1):
# residual function
resid = lambda p, x, y : sumexp(p, x) - y
# vector for initial guess
p0 = np.zeros(2 * n)
# Set new values (last two) / Use old values p[:] to refine guess
p0[-2:] = 0.1
if(n > 1): p0[:-2] = p
p, cov, infodict, mesg, ier = leastsq(resid, p0, args=(x,y), ftol=tol,
gtol=tol, maxfev=100000,
full_output=1)
residfvec = np.dot(infodict['fvec'], infodict['fvec'])
if(residfvec < tol):
print('resid', residfvec)
fail = False
break
if(fail):
raise Exception("Unable to converge decayingFunction to the desired "
+ "tolerance%3.6e with the given number of terms!"%(tol)
+ " Try increasing maxnterms or decreasing tol! "
+ "Reached%3.6e"%(residfvec))
return p
[docs]def sumexp(p, x):
"""
Returns :math:`y = \sum_n a_n b_n^{r-1}`.
**Arguments**
p : 1d array
p is an array of length 2n with n pairs. The first entry of
each pair is the coefficient :math:`a_n`, the second entry of each
pair is the exponential decay value :math:`b_n`. For example
:math:`a_0 = p_0, b_0 = p_1`.
x : 1d array
contains the points x to evaluate the function f(x). Array is of
some length l.
"""
val = np.double(0.0)
nterms = len(p) // 2
for n in range(nterms):
val = val + p[2 * n] * (p[2 * n + 1]**(x - 1))
return val
[docs]def hash_dict(dic):
"""
Generate the hash of a dictionary. Since the ordering of the entries of
`str(dic)` is always different, we build pairs of keys and values and
sort them alphabetical.
**Arguments**
dic : Dictionary
get hash of this dictionary
**Details**
Nested structures might still lead to errors in this function. Dictionaries
inside `dic` are solved by recursion calling `hash_dict`. List of
dictionaries or similar structures do not consider the dictionaries inside
the list properly.
"""
# Build two lists with keys and values
keys = []
values = []
for key in dic:
keys.append(str(key))
if(isinstance(dic[key], dict)):
values.append(hash_dict(dic[key]))
elif(hasattr(dic[key], 'get_hash')):
values.append(dic[key].get_hash())
else:
values.append(str(dic[key]))
# Sort the entries and build the string
perm = sorted(range(len(keys)), key=keys.__getitem__)
hstr = ""
for ii in perm:
hstr += keys[ii] + values[ii]
# Build the hash
hh = hashlib.sha512()
hh.update(str.encode(hstr))
return hh.hexdigest()
[docs]def bose_einstein_stat(Ek, Te, kb, mu):
"""
Return probability according to Bose-Einstein statistics.
**Arguments**
Ek : float
Energy for requested probability.
Te : float
Temperature of the system.
kb : float
Boltzmann constant.
mu : float
Chemical potential of the system.
"""
if(Te == 0.0):
return 0.0
return 1.0 / (np.exp((Ek - mu) / kb / Te) - 1.0)
[docs]def fermi_dirac_stat(Ek, Te, kb, mu):
"""
Return probability according to Bose-Einstein statistics.
**Arguments**
Ek : float
Energy for requested probability.
Te : float
Temperature of the system.
kb : float
Boltzmann constant.
mu : float
Chemical potential of the system.
"""
if(Te == 0.0):
return 0.0
return 1.0 / (np.exp((Ek - mu) / kb / Te) + 1.0)
[docs]class _MPOTerm():
"""
Abstract class providing general class methods for MPO terms.
The following class methods have to be provided for each new MPO term to work
seemlessly inside of the whole rule set for MPS calculations:
* ``__init__``
* ``AddTerm``
* ``write``
* ``print_term`` (if it should appear in the printed equation of the
Hamiltonian)
For calculations with the ED module inside openMPS the following class methods
are necessary:
* ``build_nosymm``
* ``build_symm``
"""
[docs] def check_ti(self, params):
"""
Check if term is translational invariance.
For translational invariance we check if the coupling have the
attribute ``__len__``.
**Arguments**
params : dictionary
The dictionary setting up a simulation. It coutains the couplings
which are checked here.
"""
for key in self.hp:
if(key != 'one'):
# Check variational invariance
if(hasattr(params[key], '__len__')):
return False
return True
[docs] def __len__(self):
"""
Return number of terms from hamiltonian parameters.
"""
return len(self.hp)
[docs] def __getitem__(self, key):
"""
Return the class variable according to a string
**Arguments**
key : str
Identifier string for class variable in MPO term.
"""
return self.__dict__[key]
[docs] def get_hash(self):
"""
Return a hash of the object.
"""
hstr = ''
for key in sorted(self.__dict__.keys()):
if(key in ['Operators', 'etmpo']): continue
hstr += key + str(self.__dict__[key])
hh = hashlib.sha512()
hh.update(str.encode(hstr))
return hh.hexdigest()
[docs] def ignored(self, optargs, kwargs):
"""
Print warnings about ignored optional keyword parameters:
**Arguments**
optargs : list of strings
contains all the string identifiers of the keyword arguments
considered.
kwargs : dictionary
contains all the keywords arguments
"""
for key in kwargs:
if(key not in optargs):
print("The optional argument " + key + " is not considered.")
[docs] def required(self, args, kwargs):
"""
Raise exception for required arguments not present in the optional
arguments.
**Arguments**
optargs : list of strings
contains all the string identifiers of the keyword arguments
required.
kwargs : dictionary
contains all the keywords arguments
"""
for elem in args:
if(elem not in kwargs):
raise Exception("Required argument not passed to AddMPOTerm.")
[docs] def default(self, defargs, kwargs):
"""
Set the default arguments for a term if not passed as optional
arguments in ``AddTerm``.
defargs : dict
keywords and default value for default arguments. Only
set if not present in ``kwargs``
kwargs : dict
dictionary of keywords arguments as passed to ``AddTerm``
"""
for key in defargs:
if(key not in kwargs):
kwargs[key] = defargs[key]
[docs] def herm_conj_term(self, term0, term1):
"""
Check on the hermitian conjugated term and return parameters if
existant as left operator, right operator, left weight * right weight.
**Arguments**
term0 : str
string identifier for left term
term1 : str
string identifier for right term
"""
# Check each operator for hermiticity
lh = (np.conj(self.Operators[term0].T) == self.Operators[term0]).all()
rh = (np.conj(self.Operators[term1].T) == self.Operators[term1]).all()
if(lh and rh): return None, None, None
[Lt, wlt] = self.Operators.check(np.conj(self.Operators[term0].T))
[Rt, wrt] = self.Operators.check(np.conj(self.Operators[term1].T))
return Lt, Rt, wlt * wrt
[docs] def phased_operator(self, term0, callingfunc):
"""
Generate the phased right operator and return name and weight
**Arguments**
term0 : str
string identifier for the left operator where the phase is applied
from the right.
callingfunc : str
name of the calling function to display error message.
"""
if('FermiPhase' not in self.Operators):
raise MPSOperatorNotFound(('FermiPhase', callingfunc))
PhasedR = np.dot(self.Operators[term0], self.Operators['FermiPhase'])
[name, weight] = self.Operators.check(PhasedR)
return name, weight
[docs] def phased_herm_conj_term(self, PhasedR, term1):
"""
Generate the hermitian conjugated terms for a phased pair of operators.
Return left operator, right operator, combined weight
**Arguments**
PhasedR : str
string identifier for phased operator acting on the left site.
term1 : str
string identifier for the operator acting on the right site.
"""
[leftname, lw] = self.Operators.check(np.conj(self.Operators[PhasedR].T))
[oname, ow] = self.Operators.check(np.conj(self.Operators[term1].T))
return leftname, oname, lw * ow
[docs] def decompose(self, terms, phase, gen2):
"""
Run checks and iterate over decomposition of operators and their
weight. Returns tuple (left Operators, right Operator, left weight,
right weight).
**Arguments**
term : str or 2d numpy array
Matrix acting on two sites, either in part of OperatorList or
directly numpy 2d matrix.
Phase : bool
Flag for Fermi phase.
**Details**
Checks contain hermiticity of operator and no phase term.
"""
if(phase is not None):
raise Exception("Cannot use Phase together with "
+ "decomposition request!")
if(isinstance(terms, str)):
if(not (terms in self.Operators)):
raise MPSOperatorNotFound((terms, "AddMPOTerm"))
Op2 = self.Operators[terms]
else:
Op2 = terms
if(not (np.conj(Op2.T) == Op2).all()):
raise Exception("Operator passed in to decomposition request"
+ " is not Hermitian!")
OperatorList = DecomposeOperator(Op2, gen2)
for elem in OperatorList:
[Ol, aw1] = self.Operators.check(elem[0])
[Or, aw2] = self.Operators.check(elem[1])
yield Ol, Or, aw1, aw2
[docs] def hparam(self, ii, param, quenchtime):
"""
Return actual parameter as float for a term in the
Hamiltonian multiplied with the weight.
**Arguments**
ii : int
index of the Hamiltonian parameter.
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).
"""
hpstr = self.hp[ii]
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 * self.weight[ii]
[docs] def get_sparse_dense_func(self, sparse, quenchtime):
"""
Set functions depending on sparse or dense structure of algorithms.
The functions returned are the kronecker product, the identity, and
the getting the operator.
**Operators**
sparse : bool
If true, sparse methods from ``scipy.sparse`` are used,
otherwise dense methods from ``numpy``.
quenchtime : ``None`` or tuple
For ``None`` sparse matrices are in ``csc`` format, otherwise
in ``csr`` format. No effect on dense matrices.
"""
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
if(fmt == 'csc'):
get_op = lambda op : sp.csc_matrix(self.Operators[op])
else:
get_op = lambda op : sp.csr_matrix(self.Operators[op])
else:
kron = np.kron
eye = np.eye
get_op = lambda op : self.Operators[op]
return kron, eye, get_op
[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]class SiteTerm(_MPOTerm):
"""
MPO rule set for site terms of the type
:math:`weight \cdot hparam \cdot \sum_i Op_i`.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
the operators on the local Hilbert space used.
**Variables**
Op : list
The ii-th entry contains the string identifiers for the
operators of the ii-th site term.
weight : list
The ii-th float in the list contains the weight for the
ii-th site term.
hp : list
The ii-th string in the list contains links to the Hamiltonian
paramater with the coupling used for the ii-th site term.
**Details**
To add a SiteTerm, you have to pass the following arguments to
:py:func:`mpo.MPO.AddMPOTerm`:
termtype : str
String must be equal to ``site``.
terms : str
string is the name of the operator acting locally on all sites.
weight : float, optional
Overall weight for the rule
Default to 1.0.
hparam : str, optional
String identifying the Hamiltonian parameters allowing for quenching
the parameter or changing it for different static simulations.
Default to 1.0 (``one``).
"""
def __init__(self, Operators):
self.Operators = Operators
self.Op = []
self.weight = []
self.hp = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add a new SiteTerm. For arguments see class description
:py:class:`mpoterm.SiteTerm`.
"""
# Print warnings / check required arguments / set default values
self.ignored(['weight', 'hparam'], kwargs)
self.default({'weight' : 1.0, 'hparam' : 'one'}, kwargs)
if(terms not in self.Operators):
raise MPSOperatorNotFound((terms, "SiteTerm.AddTerm"))
if((np.conj(self.Operators[terms].T) != self.Operators[terms]).any()):
raise Exception("Operator passed into site term is not "
+ "Hermitian!")
# Add site term
self.Op.append(terms)
self.weight.append(kwargs['weight'])
self.hp.append(kwargs['hparam'])
return
[docs] def write(self, fh, IntHparam, IntOperators, Params):
"""
Write the information for the fortran interface.
**Arguments**
fh : filehandle
Open filehandle of file to write to.
IntHparam : dictionary
mapping the string identifiers for Hamiltonian parameters to
integer identifiers used within fortran.
IntOperators : dictionary
mapping the string identifiers for the operators to integer
identifiers used within fortran.
Params : dictionary
containing the parameters for one simulation
"""
fh.write('%16i\n'%(len(self)))
for ii in range(len(self)):
hpw = 0 if(self.hp[ii] == 'one') else IntHparam[self.hp[ii]]
Hstring = '%30.15E'%(self.weight[ii])
Hstring += '%16i'%(hpw)
Hstring += '%16i'%(IntOperators[self.Op[ii]])
fh.write(Hstring + '\n')
return
[docs] def print_term(self):
"""
Print equation for site terms.
"""
if(len(self) == 0): return ''
tstring = ''
for ii in range(len(self)):
tstring += str(self.weight[ii]) + ' * ' + self.hp[ii] + ' * '
tstring += '\sum_i ' + self.Op[ii] + '_i'
if(ii + 1 != len(self)): tstring += ' + '
return tstring + '\n'
[docs] def build_nosymm(self, param, sparse, quenchtime, ll, ld, *args):
"""
Build the Hamiltonian for a system without symmetries. This iterator
return the matrix on the complete Hilbert space.
**Arguments**
param : dictionary
dictionary with the setup for one simulation.
sparse : bool
Use sparse methods to build the whole Hamiltonian.
quenchtime : ``None`` or tuple of integer, float
If ``None``, access static coupling. If tuple, first entry
contains the index of the quench, the second entry contains
the time to evaluate the coupling.
ll : int
System size
ld : int
Local dimension for one site.
pbc : bool
used to identify periodic boundary conditions. If not referenced
for the terms, \*args is used.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Op = get_op(self.Op[ii])
for jj in range(ll):
yield kron(kron(eye(ld**jj), coupl[jj] * Op),
eye(ld**(ll - 1 - jj)))
[docs] def build_symm(self, param, SymmSec, quenchtime, ll, *args):
"""
Build the Hamiltonian for a system with symmetries. This iterator
returns a tuple of the row index, column index, and value of the
matrix element.
**Arguments**
param : dictionary
dictionary with the setup for one simulation.
SymmSec : instance of :py:class:`exactdiag.SymmetrySector`
Allows the mapping to the basis states according to the
symmetry.
quenchtime : ``None`` or tuple of integer, float
If ``None``, access static coupling. If tuple, first entry
contains the index of the quench, the second entry contains
the time to evaluate the coupling.
ll : int
System size
pbc : bool
used to identify periodic boundary conditions. If not referenced
for the terms, \*args is used.
"""
get_op = lambda opstr : self.Operators[opstr]
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Op = self.Operators[self.Op[ii]]
colidx = self._col_ind_dic(Op)
for jj in range(ll):
for k1 in range(len(SymmSec)):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
skip = False
try:
Sel = colidx[idx1[jj]]
except KeyError:
skip = True
if(skip): continue
for ijk in Sel:
idx2[jj] = ijk[0]
elem = coupl[jj] * ijk[1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
yield (k1, k2, elem)
[docs]class BondTerm(_MPOTerm):
"""
Bond term :math:`weight \cdot hparam \cdot \sum_i Op[1]_i Op[2]_{i+1}`
is added to the to the MPO. If Phase is present then
:math:`\sum_i [ (Op[1]_i P) Op[2]_{i+1}+(P Op[1]_i^{\dagger})`
:math:`Op[2]_{i+1}^{\dagger}]`
is added. This routine will also add the Hermitian conjugate unless
both operators are their own hermitian conjugates. If only one operator
is passed then a decomposition request is assumed.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
the operators on the local Hilbert space used.
**Variables**
Op : list
The ii-th entry contains the two string identifiers for the
left and right operators of the ii-th bond term.
weight : list
The ii-th float in the list contains the weight for the
ii-th bond term.
hp : list
The ii-th string in the list contains links to the Hamiltonian
paramater with the coupling used for the ii-th bond term.
**Details**
To add a BondTerm, you have to pass the following arguments to
:py:func:`mpo.MPO.AddMPOTerm`:
termtype : str
String must be equal to ``bond``.
terms : list of str
strings are the name of the operators acting on all neighboring sites.
weight : float (optional keyword argument)
Overall weight for the rule
Default to 1.0.
hparam : str (optional keyword argument)
String identifying the Hamiltonian parameters allowing for quenching
the parameter or changing it for different static simulations.
Default to 1.0 (``one``).
Phase : boolean (optional keyword argument)
Add a Fermi phase operator originating from the Jordan-Wigner
transformation, e.g. for fermionic systems. Due to the Jordan-Wigner
transformation, the order of the operator changes the Hamiltonian.
The default order is of type ``['fdagger', 'f'].
Default to ``False``.
gen2 : list of 1d np array or ``None``
The diagonal entries of the diagonal symmetry generator of a two
site operator used when decomposition request is made. Not
referenced if two operators are passed to the bond term.
"""
def __init__(self, Operators):
self.Operators = Operators
self.Op = []
self.weight = []
self.hp = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add a new BondTerm. For arguments see class description
:py:class:`mpoterm.BondTerm`.
"""
# Print warnings / check required arguments / set default values
self.ignored(['weight', 'hparam', 'Phase', 'addcc', 'gen2'], kwargs)
self.default({'weight' : 1.0, 'hparam' : 'one', 'Phase' : None,
'addcc': True, 'gen2' : None}, kwargs)
has_phase = False if(kwargs['Phase'] is None) else kwargs['Phase']
if(not (isinstance(terms, list) or isinstance(terms, tuple))):
# Decomposition request
# .....................
for Ol, Or, wl, wr in self.decompose(terms, kwargs['Phase'],
kwargs['gen2']):
self.add([Ol, Or], kwargs['weight'] * wl * wr, kwargs['hparam'])
return
# Local operators as product
# --------------------------
if(terms[0] not in self.Operators):
raise MPSOperatorNotFound((terms[0], "BondTerm.AddTerm"))
if(terms[1] not in self.Operators):
raise MPSOperatorNotFound((terms[1], "BondTerm.AddTerm"))
if(not has_phase):
# No phase operator
# .................
self.add(terms, kwargs['weight'], kwargs['hparam'])
# Check if hermitian conjugate should be added
Lt, Rt, wlr = self.herm_conj_term(terms[0], terms[1])
if(kwargs['addcc'] and (Lt is not None)):
self.add([Lt, Rt], kwargs['weight'] * wlr, kwargs['hparam'])
return
# Phase operator
# ..............
if(('dagger' in terms[1]) and ('dagger' not in terms[0])):
print("Warning: Ordering of fermionic operators matters for the " +
"definition of phased operators. The default order is of " +
"the type ``['fdagger', 'f']``. It looks this order is not " +
"considered here. Make sure you intended to use the " +
"present order.")
rightname, rightwei = self.phased_operator(terms[0], 'BondTerm.AddTerm')
self.add([rightname, terms[1]], rightwei * kwargs['weight'],
kwargs['hparam'])
leftname, oname, low = self.phased_herm_conj_term(rightname, terms[1])
low *= rightwei * kwargs['weight']
self.add([leftname, oname], low, kwargs['hparam'])
return
[docs] def add(self, terms, weight, hparam):
"""
Adding BondTerm for internal use.
**Arguments**
terms: list of two string
Two string identifiers for the operators as list.
weight: float
Weight for the whole bond term.
hparam: str
String identifier for the Hamiltonian parameter.
"""
self.Op.append(terms)
self.weight.append(weight)
self.hp.append(hparam)
return
[docs] def write(self, fh, IntHparam, IntOperators, Params):
"""
Write the information for fortran interfaces. For description of
arguments see :py:func:`mpoterm.SiteTerm.write`.
"""
fh.write('%16i\n'%(len(self)))
for ii in range(len(self)):
hpw = 0 if(self.hp[ii] == 'one') else IntHparam[self.hp[ii]]
Hstring = '%30.15E'%(self.weight[ii])
Hstring += '%16i'%(hpw)
for jj in range(len(self.Op[ii])):
Hstring += '%16i'%(IntOperators[self.Op[ii][jj]])
fh.write(Hstring + '\n')
return
[docs] def print_term(self):
"""
Print equation for bond terms.
"""
if(len(self) == 0): return ''
tstring = ''
for ii in range(len(self)):
tstring += str(self.weight[ii]) + ' * ' + self.hp[ii] + ' * \sum_i'
tstring += self.Op[ii][0] + '_i ' + self.Op[ii][1] + '_{i+1}'
if(ii + 1 != len(self)): tstring += ' + '
return tstring + '\n'
[docs] def build_nosymm(self, param, sparse, quenchtime, ll, ld, pbc):
"""
Build the Hamiltonian for a system without symmetries. This iterator
return the matrix on the complete Hilbert space. For description of
arguments look into :py:func:`mpoterm.SiteTerm.build_nosymm`.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Op = [get_op(self.Op[ii][0]), get_op(self.Op[ii][1])]
for jj in range(ll - 1):
yield kron(kron(eye(ld**jj), coupl[jj] * Op[0]),
kron(Op[1], eye(ld**(ll - 2 - jj))))
if(pbc):
yield kron(Op[1], kron(eye(ld**(ll - 2)), coupl[-1] * Op[0]))
[docs] def build_symm(self, param, SymmSec, quenchtime, ll, pbc):
"""
Build the Hamiltonian for a system with symmetries. This iterator
returns a tuple of the row index, column index, and value of the
matrix element. For complete description look into
:py:func:`mpoterm.SiteTerm.build_symm`.
"""
get_op = lambda opstr : self.Operators[opstr]
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Op = [get_op(self.Op[ii][0]), get_op(self.Op[ii][1])]
# Get sparse representation
colidxl = self._col_ind_dic(Op[0])
colidxr = self._col_ind_dic(Op[1])
for jj in range(ll - 1):
for k1 in range(len(SymmSec)):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
# Build list with col indices
Sel = []
skip = False
try:
Sel.append(colidxl[idx1[jj]])
Sel.append(colidxr[idx1[jj + 1]])
except KeyError:
skip = True
if(skip): continue
for ijk in Cartesian(*tuple(Sel)):
elem = coupl[jj]
# Left operator
idx2[jj] = ijk[0][0]
elem *= ijk[0][1]
# Right operator
idx2[jj + 1] = ijk[1][0]
elem *= ijk[1][1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
yield (k1, k2, elem)
if(pbc):
for k1 in range(len(SymmSec)):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
# Build list with col indices
Sel = []
skip = False
try:
Sel.append(colidxl[idx1[ll - 1]])
Sel.append(colidxr[idx1[0]])
except KeyError:
skip = True
if(skip): continue
for ijk in Cartesian(*tuple(Sel)):
elem = coupl[ll - 1]
# Left operator on site ll-1
idx2[ll - 1] = ijk[0][0]
elem *= ijk[0][1]
# Right operator on site 0
idx2[0] = ijk[1][0]
elem *= ijk[1][1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
yield (k1, k2, elem)
[docs]class ExpTerm(_MPOTerm):
"""
Add an exponential term :math:`hparam \cdot \sum_{i<j} Operators[0]_i`
:math:`Operators[1]_j decayparam^{j-i-1}` to the MPO. If only one
operator is passed then a decomposition request is assumed.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
the operators on the local Hilbert space used.
**Variables**
Op : list
The ii-th entry contains the two string identifiers for the
left and right operators of the ii-th FiniteFunc term.
weight : list
The ii-th float in the list contains the weight for the
ii-th FiniteFunc term.
hp : list
The ii-th string in the list contains links to the Hamiltonian
paramater with the coupling used for the ii-th FiniteFunc term.
dp : list of floats
List of the strenght of the exponential decays of the ii-th term.
**Details**
To add a ExpTerm, you have to pass the following arguments to
:py:func:`mpo.MPO.AddMPOTerm`:
termtype : str
String must be equal to ``exponential``.
terms : list of str
strings are the name of the operators acting on all neighboring sites.
decayparam : float (keyword argument)
configuring the decay as 1 / decayparam**d where d is the distance.
weight : float (optional keyword argument)
Overall weight for the rule
Default to 1.0.
hparam : str (optional keyword argument)
String identifying the Hamiltonian parameters allowing for quenching
the parameter or changing it for different static simulations.
Default to 1.0 (``one``).
Phase : boolean (optional keyword argument)
Add a Fermi phase operator originating from the Jordan-Wigner
transformation, e.g. for fermionic systems. Due to the Jordan-Wigner
transformation, the order of the operator changes the Hamiltonian.
The default order is of type ``['fdagger', 'f'].
Default to ``False``.
"""
def __init__(self, Operators):
self.Operators = Operators
self.Op = []
self.weight = []
self.hp = []
self.dp = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add a new ExpTerm. For arguments see class description
:py:class:`mpoterm.ExpTerm`.
"""
# Print warnings / check required arguments / set default values
self.ignored(['weight', 'hparam', 'Phase', 'decayparam'], kwargs)
self.required(['decayparam'], kwargs)
self.default({'weight' : 1.0, 'hparam' : 'one', 'Phase' : None},
kwargs)
has_phase = False if(kwargs['Phase'] is None) else kwargs['Phase']
weight = kwargs['weight']
hparam = kwargs['hparam']
decayp = kwargs['decayparam']
if(not (isinstance(terms, list) or isinstance(terms, tuple))):
# Decomposition request
# .....................
for Ol, Or, wl, wr in self.decompose(terms, kwargs['Phase']):
self.add([Ol, Or, 'I'], weight * wl * wr, hparam, decayp)
return
# Local operators as product
# --------------------------
if(terms[0] not in self.Operators):
raise MPSOperatorNotFound((terms[0], "ExpTerm.AddTerm"))
if(terms[1] not in self.Operators):
raise MPSOperatorNotFound((terms[1], "ExpTerm.AddTerm"))
if(not has_phase):
# No phase operator
# .................
self.add([terms[0], terms[1], 'I'], weight, hparam, decayp)
# Check if hermitian conjugate should be added
Lt, Rt, wlr = self.herm_conj_term(terms[0], terms[1])
if(Lt is not None):
self.add([Lt, Rt, 'I'], weight, hparam, decayp)
return
# Phase operator
# ..............
if(('dagger' in terms[1]) and ('dagger' not in terms[0])):
print("Warning: Ordering of fermionic operators matters for the " +
"definition of phased operators. The default order is of " +
"the type ``['fdagger', 'f']``. It looks this order is not " +
"considered here. Make sure you intended to use the " +
"present order.")
rightname, rightwei = self.phased_operator(terms[0], 'ExpTerm.AddTerm')
self.add([rightname, terms[1], 'FermiPhase'], rightwei * weight,
hparam, decayp)
leftname, oname, low = self.phased_herm_conj_term(rightname, terms[1])
low *= rightwei * weight
self.add([leftname, oname, 'FermiPhase'], low, hparam, decayp)
return
[docs] def add(self, terms, weight, hparam, decayparam):
"""
Adding ExpTerm for internal use.
**Arguments**
terms: list of three string
Two string identifiers for the operators as list.
weight: float
Weight for the whole bond term.
hparam: str
String identifier for the Hamiltonian parameter.
decayparam : float
Decay parameter of the exponential rule.
"""
self.Op.append(terms)
self.weight.append(weight)
self.hp.append(hparam)
self.dp.append(decayparam)
[docs] def write(self, fh, IntHparam, IntOperators, Params):
"""
Write the information for fortran interfaces. For description of
arguments see :py:func:`mpoterm.SiteTerm.write`.
"""
fh.write('%16i\n'%(len(self)))
for ii in range(len(self)):
hpw = 0 if(self.hp[ii] == 'one') else IntHparam[self.hp[ii]]
dp = self.dp[ii]
Hstring = '%30.15E'%(self.weight[ii])
Hstring += '%30.15E'%(Params[dp] if(isinstance(dp, str)) else dp)
Hstring += '%16i'%(hpw)
for jj in range(len(self.Op[ii])):
Hstring += '%16i'%(IntOperators[self.Op[ii][jj]])
fh.write(Hstring + '\n')
[docs] def print_term(self):
"""
Print equation for exponential terms.
"""
if(len(self) == 0): return ''
tstring = ''
for ii in range(len(self)):
tstring += str(self.weight[ii]) + ' * ' + self.hp[ii] + ' * '
tstring += '\sum_{j<i}' + self.Op[ii][0] + '_j * ' + str(self.dp[ii])
tstring += '^{i-j-1} * ' + self.Op[ii][1] + '_i'
if(ii + 1 != len(self)): tstring += ' + '
return tstring + '\n'
[docs] def build_nosymm(self, param, sparse, quenchtime, ll, ld, *args):
"""
Build the Hamiltonian for a system without symmetries. This iterator
return the matrix on the complete Hilbert space. For description of
arguments look into :py:func:`mpoterm.SiteTerm.build_nosymm`.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Op = [get_op(self.Op[ii][0]), get_op(self.Op[ii][1]),
get_op(self.Op[ii][2])]
dp = self.dp[ii]
if(self.Op[ii][2] != 'I'):
# Terms have phase
for jj in range(ll - 1):
Tmp = kron(eye(ld**(jj)), coupl[jj] * Op[0])
for kk in range(jj + 1, ll):
dist = kk - jj
Tmp2 = deepcopy(Tmp)
for nn in range(dist - 1):
Tmp2 = kron(Tmp2, dp * Op[2])
yield kron(kron(Tmp2, Op[1]),
eye(ld**(ll - jj - dist - 1)))
else:
# Terms have no phase
for jj in range(ll - 1):
Tmp = kron(eye(ld**(jj)), coupl[jj] * Op[0])
for kk in range(jj + 1, ll):
dist = kk - jj
Tmp2 = dp**(dist - 1) * Tmp
yield kron(kron(Tmp2, eye(ld**(dist - 1))),
kron(Op[1], eye(ld**(ll - jj - dist - 1))))
[docs] def build_symm(self, param, SymmSec, quenchtime, ll, *args):
"""
Build the Hamiltonian for a system with symmetries. This iterator
returns a tuple of the row index, column index, and value of the
matrix element. For complete description look into
:py:func:`mpoterm.SiteTerm.build_symm`.
"""
get_op = lambda opstr : self.Operators[opstr]
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Op = [get_op(self.Op[ii][0]), get_op(self.Op[ii][1]),
get_op(self.Op[ii][2])]
dp = self.dp[ii]
# Sparse representation of operator
Opl = self._col_ind_dic(Op[0])
Opr = self._col_ind_dic(Op[1])
# Check if phase is not identity
phase = (self.Op[ii][2] != 'I')
if(phase):
# Last sparse representation
Opp = self._col_ind_dic(dp * Op[2])
for jj in range(ll - 1):
for hh in range(jj + 1, ll):
dist = hh - jj
#Ops = [Opl] + [Opp] * (dist - 1) + [Opr]
for k1 in range(len(SymmSec)):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
# Build list with col indices
Sel = []
skip = False
try:
Sel.append(Opl[idx1[jj]])
for kk in range(1, dist):
Sel.append(Opp[idx1[jj + kk]])
Sel.append(Opr[idx1[hh]])
except KeyError:
skip = True
if(skip): continue
for ijk in Cartesian(*tuple(Sel)):
elem = coupl[jj]
for kk in range(dist + 1):
idx2[jj + kk] = ijk[kk][0]
elem *= ijk[kk][1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
yield (k1, k2, elem)
else:
# No phase
for jj in range(ll - 1):
for hh in range(jj + 1, ll):
dist = hh - jj
for k1 in range(len(SymmSec)):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
# Build list with col indices
Sel = []
skip = False
try:
Sel.append(Opl[idx1[jj]])
Sel.append(Opr[idx1[hh]])
except KeyError:
skip = True
if(skip): continue
for ijk in Cartesian(*tuple(Sel)):
elem = coupl[jj] * dp**(dist - 1)
# Left operator
idx2[jj] = ijk[0][0]
elem *= ijk[0][1]
# Right operator
idx2[hh] = ijk[1][0]
elem *= ijk[1][1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
yield (k1, k2, elem)
[docs]class FiniteFunc(_MPOTerm):
"""
Add a FiniteFunction term
:math:`weight \cdot hparam \cdot \sum_{i,j<=r_c} f[i-j] Op[1]_i Op[2]_{j}`
to the MPO. If only one operator is passed then a decomposition request
is assumed.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
the operators on the local Hilbert space used.
**Variables**
Op : list
The ii-th entry contains the two string identifiers for the
left and right operators of the ii-th FiniteFunc term.
weight : list
The ii-th float in the list contains the weight for the
ii-th FiniteFunc term.
hp : list
The ii-th string in the list contains links to the Hamiltonian
paramater with the coupling used for the ii-th FiniteFunc term.
f : list of 1d vectors
List of coupling up to the specified range.
r_c : list of ints
Range of the ii-th FiniteFunc term.
Phase : list of strings
List of the phase operator inbetween the left and right operator.
This is the identity ``'I'`` if no phase present.
**Details**
To add a FiniteFunc term, you have to pass the following arguments to
:py:func:`mpo.MPO.AddMPOTerm`
termtype : str
String must be equal to ``FiniteFunc``.
terms : str
string is the name of the operator acting locally on all sites.
f : 1d array (keyword argument)
Containing the weight for a with distance for a range specified
through the length of the array.
weight : float (optional keyword argument)
Overall weight for the rule
Default to 1.0.
hparam : str (optional keyword argument)
String identifying the Hamiltonian parameters allowing for quenching
the parameter or changing it for different static simulations.
Default to 1.0 (``one``).
Phase : boolean (optional keyword argument)
Add a Fermi phase operator originating from the Jordan-Wigner
transformation, e.g. for fermionic systems. Due to the Jordan-Wigner
transformation, the order of the operator changes the Hamiltonian.
The default order is of type ``['fdagger', 'f'].
Default to ``False``.
"""
def __init__(self, Operators):
self.Operators = Operators
self.Op = []
self.weight = []
self.hp = []
self.f = []
self.r_c = []
self.Phase = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add a new FiniteFunc term. For arguments see class description
:py:class:`mpoterm.FiniteFunc`.
"""
self.ignored(['weight', 'hparam', 'Phase', 'f'], kwargs)
self.required(['f'], kwargs)
self.default({'weight' : 1.0, 'hparam' : 'one', 'Phase' : False}, kwargs)
has_phase = False if(kwargs['Phase'] is None) else kwargs['Phase']
# Convert to numpy array / set phase operator
kwargs['f'] = np.array(kwargs['f'])
PhaseOp = 'FermiPhase' if(kwargs['Phase']) else 'I'
if(not (isinstance(terms, list) or isinstance(terms, tuple))):
# Decomposition request
# .....................
for Ol, Or, wl, wr in self.decompose(terms, kwargs['Phase']):
self.add(terms, kwargs['weight'], kwargs['hparam'],
kwargs['f'], PhaseOp)
return
# Local operators as product
# --------------------------
if(not (terms[0] in self.Operators)):
raise MPSOperatorNotFound((terms[0], "FiniteFunc.AddTerm"))
if(not (terms[1] in self.Operators)):
raise MPSOperatorNotFound((terms[1], "FiniteFunc.AddTerm"))
if(not has_phase):
# No phase operator
# .................
self.add(terms, kwargs['weight'], kwargs['hparam'], kwargs['f'],
PhaseOp)
# Check if hermitian conjugate should be added
Lt, Rt, wlr = self.herm_conj_term(terms[0], terms[1])
if(Lt is not None):
self.add([Lt, Rt], kwargs['weight'] * wlr, kwargs['hparam'],
kwargs['f'], PhaseOp)
return
# Phase operator
# ..............
if(('dagger' in terms[1]) and ('dagger' not in terms[0])):
print("Warning: Ordering of fermionic operators matters for the " +
"definition of phased operators. The default order is of " +
"the type ``['fdagger', 'f']``. It looks this order is not " +
"considered here. Make sure you intended to use the " +
"present order.")
rightname, rightwei = self.phased_operator(terms[0], 'FiniteFunc.AddTerm')
self.add([rightname, terms[1]], rightwei * kwargs['weight'],
kwargs['hparam'], kwargs['f'], PhaseOp)
leftname, oname, low = self.phased_herm_conj_term(rightname, terms[1])
low *= rightwei * kwargs['weight']
self.add([leftname, oname], low, kwargs['hparam'], kwargs['f'], PhaseOp)
return
[docs] def add(self, terms, weight, hparam, farray, Phase):
"""
Adding FiniteFunc term for internal use.
**Arguments**
terms: list of two string
Two string identifiers for the operators as list.
weight: float
Weight for the whole bond term.
hparam: str
String identifier for the Hamiltonian parameter.
farray : 1d vector
Coupling as a function of distance. Specifies the
range through the length of the array.
Phase : str
String identifier for the phase operator.
"""
self.Op.append(terms)
self.weight.append(weight)
self.hp.append(hparam)
self.f.append(farray)
self.r_c.append(farray.shape[0])
self.Phase.append(Phase)
[docs] def write(self, fh, IntHparam, IntOperators, Params):
"""
Write the information for fortran interfaces. For description of
arguments see :py:func:`mpoterm.SiteTerm.write`.
"""
fh.write('%16i\n'%(len(self)))
for ii in range(len(self)):
hpw = 0 if(self.hp[ii] == 'one') else IntHparam[self.hp[ii]]
Hstring = '%30.15E'%(self.weight[ii])
Hstring += '%16i'%(hpw)
Hstring += '%16i'%(self.r_c[ii])
Hstring += '%16i'%(IntOperators[self.Op[ii][0]])
Hstring += '%16i'%(IntOperators[self.Op[ii][1]])
Hstring += '%16i'%(IntOperators[self.Phase[ii]])
fh.write(Hstring + '\n')
Hstring = ''
for jj in range(self.r_c[ii]):
Hstring += '%30.15E'%(self.f[ii][jj])
fh.write(Hstring + '\n')
[docs] def print_term(self):
"""
Print equation for FiniteFunction terms.
"""
if(len(self) == 0): return ''
tstring = ''
for ii in range(len(self)):
tstring += str(self.weight[ii]) + ' * ' + self.hp[ii] + ' * '
tstring += '\sum_{i,i+1<=j<=i+' + str(self.r_c[ii]) + '} f^{'
tstring += str(ii) + '}_{j-i}' + self.Op[ii][0] + '_i * '
tstring += self.Op[ii][1] + '_j'
if(ii + 1 != len(self)): tstring += ' + '
tstring += '\nWhere\n'
for ii in range(len(self)):
tstring += 'f^{' + str(ii) + '}= [ '
for jj in self.f[ii]:
tstring += str(jj) + ' '
tstring += ']\n'
return tstring
[docs]class InfiniteFunc(_MPOTerm):
"""
Add a term to the Hamiltonian
:math:`weight \cdot hparam \cdot \sum_{i<j} f(i-j) Operators[0]_i`
:math:`Operators[1]_j` where :math:`f(x)` is specified as a function
handle on input. This is done by finding the series of weighted
exponentials :math:`\sum_{a} c_a e^{-b x}`
which best approximates this function and then adding those terms to
the exponential part of the MPO. If only one operator is passed then a
decomposition request is assumed.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
the operators on the local Hilbert space used.
ExpTermMPO : instance of :py:class:`mpoterm.ExpTerm`
Add the fitted exponential terms to this instance of a class.
PostProcess : bool
If ``True`` the functions are not fitted. If ``False`` fitting
procedure expresses InfiniteFunc as ExpTerms.
**Variables**
etmpo : instance of :py:class:`mpoterm.ExpTerm`
Add the fitted exponential terms to this instance of a class. Passed
to ``__init__`` as ``ExpTermMPO``.
Op : list
The ii-th entry contains the two string identifiers for the
left and right operators of the ii-th bond term.
weight : list
The ii-th float in the list contains the weight for the
ii-th bond term.
hp : list
The ii-th string in the list contains links to the Hamiltonian
paramater with the coupling used for the ii-th bond term.
yvals : list of vectors
The ii-th vector specifies the coupling with distance.
**Details**
To add a InfiniteFunc term, you have to pass the following arguments to
:py:func:`mpo.MPO.AddMPOTerm`:
termtype : str
String must be equal to ``InfiniteFunction``.
terms : list of str
strings are the name of the operators acting on all neighboring sites.
weight : float (optional keyword argument)
Overall weight for the rule
Default to 1.0.
hparam : str (optional keyword argument)
String identifying the Hamiltonian parameters allowing for quenching
the parameter or changing it for different static simulations.
Default to 1.0 (``one``).
Phase : boolean (optional keyword argument)
Add a Fermi phase operator originating from the Jordan-Wigner
transformation, e.g. for fermionic systems. Due to the Jordan-Wigner
transformation, the order of the operator changes the Hamiltonian.
The default order is of type ``['fdagger', 'f'].
Default to ``False``.
maxnterms : int (optional keyword argument)
Number of exponentials maximally allowed to fit InfiniteFunction
Default to 16.
tol : float (optional keyword argument)
This sets the tolerance for the fit of the infinite function
with the exponentials.
Default to 1e-6
func : callable (optional keyword argument)
Function returning the coupling as a function for the distance.
Combined with the argument ``L``, it may be used to specify the
InfiniteFunction
Default to ``None``
L : int (optional keyword argument)
range for fitting the exponentials to the function ``func``
Default to ``None``
x : 1d vector (optional keyword argument)
The values of the function to be fitted can be given as two vector
x and y = f(x). These are the x-values.
Default to ``None``
y : 1d vector (optional keyword argument)
The values of the function to be fitted can be given as two vector
x and y = f(x). These are values of the function.
Default to ``None``
ED : bool (optional keyword argument)
It can only be used to indicate for the EDLib not to fit an
InfiniteFunction with exponentials.
Default to ``False``
Although func, L, x, y are all marked as optional, either x and y or
func and L have to be provided.
"""
def __init__(self, Operators, ExpTermMPO, PostProcess):
self.Operators = Operators
self.etmpo = ExpTermMPO
self.PostProcess = PostProcess
self.Op = []
self.weight = []
self.hp = []
self.yvals = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add a new InfiniteFunc. For arguments see class description
:py:class:`mpoterm.InfiniteFunc`.
"""
# Print warnings / check required arguments / set default values
args = ['weight', 'hparam', 'Phase', 'maxnterms', 'tol', 'func', 'L',
'x', 'y', 'ED']
defaults = {'weight' : 1.0, 'hparam' : 'one', 'Phase' : None,
'maxnterms' : 16, 'tol' : 1e-6, 'func' : None,
'L' : None, 'x' : None, 'y' : None, 'ED' : False}
self.ignored(args, kwargs)
self.default(defaults, kwargs)
if(self.PostProcess):
# Need to add dummy term for checking hparam in AddQuench
self.etmpo.AddTerm(terms, hparam=kwargs['hparam'], decayparam=1.0,
weight=0.0, Phase=kwargs['Phase'])
return
# Check x, y, func, and L directly
x = kwargs['x']
y = kwargs['y']
L = kwargs['L']
func = kwargs['func']
if((x is None) and (y is None) and (func is not None)
and (L is not None)):
# Generate values from function
x = np.linspace(1, L, L)
y = func(x)
elif((x is not None) and (y is not None)):
pass
else:
raise Exception("Either func and L or x and y should be supplied to"
+ " AddInfiniteFunctionMPOterm!")
if(kwargs['ED']):
self.AddTermED(terms, kwargs['weight'], kwargs['hparam'],
kwargs['Phase'], x, y)
return
# MPS InfiniteFunc are fitted by series of exponentials
p = fittoexpSum(x, y, kwargs['maxnterms'], kwargs['tol'])
nn = len(p) // 2
print('Number of exponentials used to fit infinite function', nn)
fitparams = []
for ii in range(nn):
fitparams.append([p[2 * ii], p[2 * ii + 1]])
# Local operators or decomposition request will be treated inside
# of AddexponentialMPOterm
kwargs2 = {'hparam' : kwargs['hparam'], 'Phase' : kwargs['Phase']}
for kk in fitparams:
kwargs2['weight'] = kwargs['weight'] * kk[0]
kwargs2['decayparam'] = kk[1]
self.etmpo.AddTerm(terms, **kwargs2)
return
[docs] def AddTermED(self, terms, weight, hparam, Phase, xval, yval):
"""
For the exact diagonlization code we can use the infinite functions
directly.
**Arguments**
terms: list of two string
Two string identifiers for the operators as list.
weight: float
Weight for the whole bond term.
hparam: str
String identifier for the Hamiltonian parameter.
Phase : str
String identifier for the phase operator.
xval : 1d vector
values with the x-coordinates for the yval = f(xvals). Has to
start at 1 and have an increment 1.
yval : 1d vector
Values with the distance dependent couplting.
"""
if((xval[0] != 1) or np.any(xval[1:] - xval[:-1] != 1)):
raise Exception('x values need to start at 1 and have state 1.')
has_phase = False if(Phase is None) else Phase
if(not (isinstance(terms, list) or isinstance(terms, tuple))):
# Decomposition request
# .....................
for Ol, Or, wl, wr in self.decompose(terms, Phase):
self.add([Ol, Or, 'I'], weight * wl * wr, hparam, yval)
return
# Local operators as product
if(terms[0] not in self.Operators):
raise MPSOperatorNotFound((terms[0], "InfiniteFunc.AddTermED"))
if(terms[1] not in self.Operators):
raise MPSOperatorNotFound((terms[1], "InfiniteFunc.AddTermED"))
if(not has_phase):
# No phase operator
# .................
self.add([terms[0], terms[1], 'I'], weight, hparam, yval)
# Check if hermitian conjugate should be added
Lt, Rt, wlr = self.herm_conj_term(terms[0], terms[1])
if(Lt is not None):
self.add([Lt, Rt, 'I'], weight, hparam, yval)
return
# Phase operator
# ..............
if(('dagger' in terms[1]) and ('dagger' not in terms[0])):
print("Warning: Ordering of fermionic operators matters for the " +
"definition of phased operators. The default order is of " +
"the type ``['fdagger', 'f']``. It looks this order is not " +
"considered here. Make sure you intended to use the " +
"present order.")
rightname, rightwei = self.phased_operator(terms[0], 'ExpTerm.AddTerm')
self.add([rightname, terms[1], 'FermiPhase'], rightwei * weight,
hparam, yval)
leftname, oname, low = self.phased_herm_conj_term(rightname, terms[1])
low *= rightwei * weight
self.add([leftname, oname, 'FermiPhase'], low, hparam, yval)
return
[docs] def add(self, terms, weight, hparam, yval):
"""
Add InfiniteFunc term, internal use only.
**Arguments**
terms: list of two string
Three string identifiers for the operators as list. Order is
left operator, right operator, phase operator.
weight: float
Weight for the whole bond term.
hparam: str
String identifier for the Hamiltonian parameter.
yval : 1d vector
Values with the distance dependent couplting.
"""
self.Op.append(terms)
self.weight.append(weight)
self.hp.append(hparam)
self.yvals.append(yval)
[docs] def print_term(self):
"""
Print equation for InfiniteFunction terms.
"""
if(len(self) == 0): return ''
print('InfiniteFunction terms not printed in equation yet.')
return ''
[docs] def build_nosymm(self, param, sparse, quenchtime, ll, ld, *args):
"""
Build the Hamiltonian for a system without symmetries. This iterator
return the matrix on the complete Hilbert space. For description of
arguments look into :py:func:`mpoterm.SiteTerm.build_nosymm`.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Op = [get_op(self.Op[ii][0]), get_op(self.Op[ii][1]),
get_op(self.Op[ii][2])]
yvals = self.yvals[ii]
if(self.Op[ii][2] != 'I'):
# Terms have phase
for jj in range(ll - 1):
Tmp = kron(eye(ld**jj), coupl[jj] * Op[0])
for kk in range(jj + 1, ll):
dist = kk - jj
Tmp2 = deepcopy(Tmp) * yvals[dist - 1]
for nn in range(dist - 1):
Tmp2 = kron(Tmp2, Op[2])
yield kron(kron(Tmp2, Op[1]),
eye(ld**(ll - jj - dist - 1)))
else:
# Terms have no phase
for jj in range(ll - 1):
Tmp = kron(eye(ld**jj), coupl[jj] * Op[0])
for kk in range(jj + 1, ll):
dist = kk - jj
Tmp2 = Tmp * yvals[dist - 1]
yield kron(kron(Tmp2, eye(ld**(dist - 1))),
kron(Op[1], eye(ld**(ll - jj - dist - 1))))
[docs] def build_symm(self, param, SymmSec, quenchtime, ll, *args):
"""
Build the Hamiltonian for a system with symmetries. This iterator
returns a tuple of the row index, column index, and value of the
matrix element. For complete description look into
:py:func:`mpoterm.SiteTerm.build_symm`.
"""
get_op = lambda opstr : self.Operators[opstr]
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Op = [get_op(self.Op[ii][0]), get_op(self.Op[ii][1]),
get_op(self.Op[ii][2])]
yvals = self.yvals[ii]
# Sparse representation of operator
Opl = self._col_ind_dic(Op[0])
Opr = self._col_ind_dic(Op[1])
# Check if phase is not identity
phase = (self.Op[ii][2] != 'I')
if(phase):
# Last sparse representation
Opp = self._col_ind_dic(Op[2])
for jj in range(ll - 1):
for hh in range(jj + 1, ll):
dist = hh - jj
#Ops = [Opl * yvals[dist - 1]] \
# + [Opp] * (dist - 1) + [Opr]
for k1 in range(len(SymmSec)):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
# Build list with col indices
Sel = []
skip = False
try:
Sel.append(Opl[idx1[jj]])
for kk in range(1, dist):
Sel.append(Opp[idx1[jj + kk]])
Sel.append(Opr[idx1[hh]])
except KeyError:
skip = True
if(skip): continue
for ijk in Cartesian(*tuple(Sel)):
elem = coupl[jj]
for kk in range(dist + 1):
idx2[jj + kk] = ijk[kk][0]
elem *= ijk[kk][1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
yield (k1, k2, elem)
else:
# No phase
for jj in range(ll - 1):
for hh in range(jj + 1, ll):
dist = hh - jj
for k1 in range(len(SymmSec)):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
# Build list with col indices
Sel = []
skip = False
try:
Sel.append(Opl[idx1[jj]])
Sel.append(Opr[idx1[hh]])
except KeyError:
skip = True
if(skip): continue
for ijk in Cartesian(*tuple(Sel)):
elem = coupl[jj] * yvals[dist - 1]
# Left operator
idx2[jj] = ijk[0][0]
elem *= ijk[0][1]
# Right operator
idx2[hh] = ijk[1][0]
elem *= ijk[1][1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
yield (k1, k2, elem)
[docs]class ProductTerm(_MPOTerm):
"""
Add a product term :math:`weight \cdot hparam \cdot \prod_i Op_i`
to the MPO.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
the operators on the local Hilbert space used.
**Variables**
Op : list
The ii-th entry contains the two string identifiers for the
left and right operators of the ii-th bond term.
weight : list
The ii-th float in the list contains the weight for the
ii-th bond term.
hp : list
The ii-th string in the list contains links to the Hamiltonian
paramater with the coupling used for the ii-th bond term.
**Details**
To add a ProductTerm, you have to pass the following arguments to
:py:func:`mpo.MPO.AddMPOTerm`:
termtype : str
String must be equal to ``product``.
terms : str
String identifier for the operator.
weight : float (optional keyword argument)
Overall weight for the rule
Default to 1.0.
hparam : str (optional keyword argument)
String identifying the Hamiltonian parameters allowing for quenching
the parameter or changing it for different static simulations.
Default to 1.0 (``one``).
"""
def __init__(self, Operators):
self.Operators = Operators
self.Op = []
self.weight = []
self.hp = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add a new BondTerm. For arguments see class description
:py:class:`mpoterm.ProductTerm`.
"""
# Print warnings / check required arguments / set default values
self.ignored(['weight', 'hparam'], kwargs)
self.default({'weight' : 1.0, 'hparam' : 'one'}, kwargs)
if(terms not in self.Operators):
raise MPSOperatorNotFound((terms, "ProductTerm.AddTerm"))
if((np.conj(self.Operators[terms].T) != self.Operators[terms]).any()):
raise Exception("Operator passed into product term is not "
+ "Hermitian!")
# Add product term
self.Op.append(terms)
self.weight.append(kwargs['weight'])
self.hp.append(kwargs['hparam'])
return
[docs] def write(self, fh, IntHparams, IntOperators, Params):
"""
Write the information for fortran interfaces. For description of
arguments see :py:func:`mpoterm.SiteTerm.write`.
"""
fh.write('%16i\n'%(len(self)))
for ii in range(len(self)):
hpw = 0 if(self.hp[ii] == 'one') else IntHparams[self.hp[ii]]
Hstring = '%30.15E'%(self.weight[ii])
Hstring += '%16i'%(hpw)
Hstring += '%16i'%(IntOperators[self.Op[ii]])
fh.write(Hstring + '\n')
[docs] def print_term(self):
"""
Print equation for product terms.
"""
if(len(self) == 0): return ''
tstring = ''
for ii in range(len(self)):
tstring += str(self.weight[ii]) + ' * ' + self.hp[ii] + ' * '
tstring += '\prod__i ' + self.Op[ii] + '_i'
if(ii + 1 != len(self)): tstring += ' + '
return tstring + '\n'
[docs]class MBStringTerm(_MPOTerm):
"""
Add a many-body string term
:math:`weight \cdot hparam \cdot \sum_i \prod_{j=1}^{n} Op[j]_{i+j-1}`
to the MPO, where Op[j] is the :math:`j^{th}` element of term. The number
of such operators n is determined from the length of the input list.
This routine does not check for Hermiticity.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
the operators on the local Hilbert space used.
**Variables**
Op : list
The ii-th entry contains the two string identifiers for the
left and right operators of the ii-th bond term.
weight : list
The ii-th float in the list contains the weight for the
ii-th bond term.
hp : list
The ii-th string in the list contains links to the Hamiltonian
paramater with the coupling used for the ii-th bond term.
**Details**
To add a MBStringTerm, you have to pass the following arguments to
:py:func:`mpo.MPO.AddMPOTerm`:
termtype : str
String must be equal to ``MBString``.
terms : list of str
String identifiers for the string of operators.
weight : float (optional keyword argument)
Overall weight for the rule
Default to 1.0.
hparam : str (optional keyword argument)
String identifying the Hamiltonian parameters allowing for quenching
the parameter or changing it for different static simulations.
Default to 1.0 (``one``).
"""
def __init__(self, Operators):
self.Operators = Operators
self.Op = []
self.weight = []
self.hp = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add a new MBString term. For arguments see class description
:py:class:`mpoterm.MBStringTerm`.
"""
self.ignored(['weight', 'hparam'], kwargs)
self.default({'weight' : 1.0, 'hparam' : 'one'}, kwargs)
if not hasattr(terms, '__len__'):
raise Exception("terms passed into AddMPOTerm for a MBString is "
+ "not a list of Operator tags!")
if(len(terms) < 2):
raise Exception("terms passed into AddMPOTerm for a MBString "
+ "should be a list with at least two elements!")
for tm in terms:
if not tm in self.Operators:
raise MPSOperatorNotFound((tm, "MBStringTerm.AddTerm"))
self.Op.append(terms)
self.weight.append(kwargs['weight'])
self.hp.append(kwargs['hparam'])
return
[docs] def write(self, fh, IntHparam, IntOperators, Params):
"""
Write the information for fortran interfaces. For description of
arguments see :py:func:`mpoterm.SiteTerm.write`.
"""
fh.write('%16i\n'%(len(self)))
for ii in range(len(self)):
hpw = 0 if(self.hp[ii] == 'one') else IntHparam[self.hp[ii]]
Hstring = '%30.15E'%(self.weight[ii])
Hstring += '%16i'%(hpw)
Hstring += '%16i'%(len(self.Op[ii])) + '\n'
for jj in range(len(self.Op[ii])):
Hstring += '%16i'%(IntOperators[self.Op[ii][jj]])
fh.write(Hstring + '\n')
[docs] def print_term(self):
"""
Print equation for MBString terms.
"""
if(len(self) == 0): return ''
tstring = ''
for ii in range(len(self)):
tstring += str(self.weight[ii]) + ' * ' + self.hp[ii] + ' * \sum_{i}'
tstring += str(self.Op[ii][0]) + '_{i}'
for jj in range(1, len(self.Op[ii])):
tstring += ' * ' + str(self.Op[ii][jj]) + '_{i+' + str(jj) + '}'
if(ii + 1 != len(self)): tstring += ' + '
return tstring + '\n'
[docs] def build_nosymm(self, param, sparse, quenchtime, ll, ld, *args):
"""
Build the Hamiltonian for a system without symmetries. This iterator
return the matrix on the complete Hilbert space. For description of
arguments look into :py:func:`mpoterm.SiteTerm.build_nosymm`.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Op = []
for opstr in self.Op[ii]:
Op.append(get_op(opstr))
nk = len(self.Op[ii])
Tmp = Op[0]
for kk in range(1, nk):
Tmp = kron(Tmp, Op[kk])
for jj in range(ll - nk + 1):
yield kron(kron(eye(ld**jj) * coupl[jj], Tmp),
eye(ld**(ll - nk - jj)))
[docs] def build_symm(self, param, SymmSec, quenchtime, ll, *args):
"""
Build the Hamiltonian for a system with symmetries. This iterator
returns a tuple of the row index, column index, and value of the
matrix element. For complete description look into
:py:func:`mpoterm.SiteTerm.build_symm`.
"""
get_op = lambda opstr : self.Operators[opstr]
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Ops = {}
mapping = []
kk = 0
jj = -1
for opstr in self.Op[ii]:
jj += 1
if(opstr != 'I'):
Ops[kk] = self._col_ind_dic(self.Operators[opstr])
mapping.append(jj)
kk += 1
nk = len(self.Op[ii])
nks = len(mapping)
for jj in range(ll - nk + 1):
for k1 in range(len(SymmSec)):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
Sel = []
skip = False
for kk in range(nks):
try:
Sel.append(Ops[kk][idx1[jj + mapping[kk]]])
except KeyError:
skip = True
break
if(skip): continue
for ijk in Cartesian(*tuple(Sel)):
elem = coupl[jj]
for kk in range(nks):
idx2[jj + mapping[kk]] = ijk[kk][0]
elem *= ijk[kk][1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
yield (k1, k2, elem)
[docs]class TTerm(_MPOTerm):
"""
MPO rule for a site 1 to k-1 interaction with site k. The interaction
strength is independent of the distance.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
the operators on the local Hilbert space used.
**Variables**
Op : list
The ii-th entry contains the two string identifiers for the
left and right operators of the ii-th bond term.
kk : list of ints
The ii-th entry contains the site which sites 1 .. (kk - 1) are
interacting to for the ii-th rule.
weight : list
The ii-th float in the list contains the weight for the
ii-th bond term.
hp : list
The ii-th string in the list contains links to the Hamiltonian
paramater with the coupling used for the ii-th bond term.
**Details**
To add an TTerm you have to pass the following arguments to
:py:func:`mpo.MPO.AddMPOTerm`:
termtype : str
must be ``tterm`` to add this term.
terms : list of two strings
strings are the names of the operators acting on the sites 1 to k-1
(first string) and the site k (second string).
kk : int (keyword argument)
All sites left of site kk interact with site k. Indexing NOT in
python standard starting from 0, but from 1.
weight : float (optional keyword argument)
Overall weight for the rule
Default to 1.0.
hparam : str (optional keyword argument)
String identifying the Hamiltonian parameters allowing for quenching
the parameter or changing it for different static simulations.
Default to 1.0 (``one``).
"""
def __init__(self, Operators):
self.Operators = Operators
self.Op = []
self.kk = []
self.weight = []
self.hp = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add a new TTerm. For arguments see class description
:py:class:`mpoterm.TTerm`
"""
print('Warning: Tterm not implemented in MPS (ED partly).')
# Print warnings / check required arguments / set default values
self.ignored(['kk', 'weight', 'hparam'], kwargs)
self.required(['kk'], kwargs)
self.default({'weight' : 1.0, 'hparam' : 'one'}, kwargs)
if(not (isinstance(terms, list) or isinstance(terms, tuple))):
# Decomposition request
# .....................
for Ol, Or, wl, wr in self.decompose(terms, None):
self.add([Ol, Or], kwargs['kk'], kwargs['weight'] * wl * wr,
kwargs['hparam'])
return
# Local operators as product (never a phase operator)
# ..........................
if(terms[0] not in self.Operators):
raise MPSOperatorNotFound((terms[0], "AddMPOTerm"))
if(terms[1] not in self.Operators):
raise MPSOperatorNotFound((terms[1], "AddMPOTerm"))
# Add term
self.add(terms, kwargs['kk'], kwargs['weight'], kwargs['hparam'])
# Check if hermitian conjugate terms should be added
Lt, Rt, wlr = self.herm_conj_term(terms[0], terms[1])
if(Lt is not None):
self.add([Lt, Rt], kwargs['kk'], kwargs['weight'] * wlr,
kwargs['hparam'])
return
[docs] def add(self, terms, kk, weight, hparam):
"""
Adding TTerm for internal use.
**Arguments**
terms: list of two string
Two string identifiers for the operators as list.
kk : int
All sites left of site kk interact with site k.
weight: float
Weight for the whole bond term.
hparam: str
String identifier for the Hamiltonian parameter.
"""
self.Op.append(terms)
self.kk.append(kk)
self.weight.append(weight)
self.hp.append(hparam)
return
[docs] def write(self, fh, IntHparam, IntOperators, Params):
"""
Write the information for fortran interfaces. For description of
arguments see :py:func:`mpoterm.SiteTerm.write`.
"""
fh.write('%16i\n'%(len(self)))
for ii in range(len(self)):
hpw = 0 if(self.hp[ii] == 'one') else IntHparams[self.hp[ii]]
Hstring = '%30.15E'%(self.weight[ii])
Hstring += '%16i'%(hpw)
Hstring += '%16i'%(self.kk[ii])
for jj in range(len(self.Op[ii])):
Hstring += '%16i'%(IntOperators[self.Op[ii][jj]])
fh.write(Hstring + '\n')
[docs] def print_term(self):
"""
Print equation for Tterms.
"""
if(len(self) == 0): return ''
print('Tterms not printed in equation yet.')
return ''
[docs] def check_ti(self, *args):
"""
TTerms are by default not translational invariant since they
depend on site kk.
"""
return (len(self) == 0)
[docs] def build_nosymm(self, param, sparse, quenchtime, ll, ld, *args):
"""
Build the Hamiltonian for a system without symmetries. This iterator
returns the matrix on the complete Hilbert space. For description
of arguments look into :py:func:mpoterm.SiteTerm.build_nosymm`.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Opl = get_op(self.Op[ii][0])
Opr = get_op(self.Op[ii][1])
# Get this kk in python notation
kk = self.kk[ii] - 1
# Filling identities between Opl and Opr
nf = kk
# Identities on the right of Opr
nr = ll - kk - 1
# Temporary operator Opr x Id
OprId = kron(Opr, eye(ld**nr))
for jj in range(kk):
nf -= 1
yield kron(kron(eye(ld**jj), coupl[jj] * Opl),
kron(eye(ld**nf), OprId))
[docs] def build_symm(self, param, SymmSec, quenchtime, ll, *args):
"""
Build the Hamiltonian for a system with symmetries. This iterator
returns a tuple of the row index, column index, and value of the
matrix element. For complete description look into
:py:func:`mpoterm.SiteTerm.build_symm`. NOT IMPLEMENTED.
"""
raise NotImplementedError
[docs]class SiteLindXX(_MPOTerm):
"""
Add a Lindblad term :math:`weight \cdot hparam \cdot \sum_i Op_i`
to the MPO.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
the operators on the local Hilbert space used.
**Variables**
Op : list
The ii-th entry contains the two string identifiers for the
left and right operators of the ii-th bond term.
weight : list
The ii-th float in the list contains the weight for the
ii-th bond term.
hp : list
The ii-th string in the list contains links to the Hamiltonian
paramater with the coupling used for the ii-th bond term.
**Details**
To add a single site Lindblad channel, you have to pass the following
arguments to :py:func:`mpo.MPO.AddMPOTerm`:
term : string
string name of the operator
weight : float, optional
general weight, e.g. a minus sign
hparam : string, optional
possibly site and time dependent coupling
default to 'one'
"""
def __init__(self, Operators):
self.Operators = Operators
self.Op = []
self.weight = []
self.hp = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add a new SiteLindXX. For arguments see class description
:py:class:`mpoterm.SiteLindXX`.
"""
# Print warnings / check required arguments / set default values
self.ignored(['weight', 'hparam'], kwargs)
self.default({'weight' : 1.0, 'hparam' : 'one'}, kwargs)
if(terms not in self.Operators):
raise MPSOperatorNotFound((terms, "SiteLindXX.AddTerm"))
# Add site term
self.Op.append(terms)
self.weight.append(kwargs['weight'])
self.hp.append(kwargs['hparam'])
return
[docs] def write(self, fh, IntHparam, IntOperators, Params):
"""
Write the information for fortran interfaces. For description of
arguments see :py:func:`mpoterm.SiteTerm.write`.
"""
fh.write('%16i\n'%(len(self)))
for ii in range(len(self)):
hpw = 0 if(self.hp[ii] == 'one') else IntHparam[self.hp[ii]]
Hstring = '%30.15E'%(self.weight[ii])
Hstring += '%16i'%(hpw)
Hstring += '%16i'%(IntOperators[self.Op[ii]])
fh.write(Hstring + '\n')
[docs] def print_term(self):
"""
Print equation for SiteLindXX terms.
"""
if(len(self) == 0): return ''
print('Lindblad terms not printed in equation yet.')
return ''
[docs] def build_nosymm(self, param, sparse, quenchtime, ll, ld, *args):
"""
Build the effective Hamiltonian for a system without symmetries. This
iterator return the matrix on the complete Hilbert space. For
description of arguments look into
:py:func:`mpoterm.SiteTerm.build_nosymm`.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Op = get_op(self.Op[ii])
Op = -0.5 * Op.conj().transpose().dot(Op)
for jj in range(ll):
yield kron(kron(eye(ld**jj), coupl[jj] * Op),
eye(ld**(ll - 1 - jj)))
[docs] def lbuild_nosymm(self, param, sparse, quenchtime, ll, ld, *args):
"""
Build the Liouville operator for a system without symmetries. This
iterators returns the matrix build on the complete Liouville space.
**Arguments**
param : dictionary
dictionary with the setup for one simulation.
sparse : bool
Use sparse methods to build the whole Hamiltonian.
quenchtime : ``None`` or tuple of integer, float
If ``None``, access static coupling. If tuple, first entry
contains the index of the quench, the second entry contains
the time to evaluate the coupling.
ll : int
System size
ld : int
Local dimension for one site.
pbc : bool
used to identify periodic boundary conditions. If not referenced
for the terms, \*args is used.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
# Terms of -0.5 Ld L x I and -0.5 I x (Ld L)^T
for elem in self.build_nosymm(param, sparse, quenchtime, ll, ld):
yield kron(elem, eye(ld**ll))
yield kron(eye(ld**ll), elem.transpose())
# Terms of Li x Ld^T = Li x Li*
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Op = get_op(self.Op[ii])
for jj in range(ll):
yield kron(kron(kron(eye(ld**(jj)), coupl[jj] * Op),
kron(eye(ld**(ll - 1)), Op.conj())),
eye(ld**(ll - 1 - jj)))
[docs] def build_symm(self, param, SymmSec, quenchtime, ll, *args):
"""
Build the effective Hamiltonian for a system with symmetries. This
iterator returns a tuple of the row index, column index, and value
of the matrix element. For complete description look into
:py:func:`mpoterm.SiteTerm.build_symm`.
"""
get_op = lambda opstr : self.Operators[opstr]
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Op = self.Operators[self.Op[ii]]
Op = -0.5 * Op.conj().transpose().dot(Op)
colidx = self._col_ind_dic(Op)
for jj in range(ll):
for k1 in range(len(SymmSec)):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
skip = False
try:
Sel = colidx[idx1[jj]]
except KeyError:
skip = True
if(skip): continue
for ijk in Sel:
idx2[jj] = ijk[0]
elem = coupl[jj] * ijk[1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
yield (k1, k2, elem)
[docs] def lbuild_symm(self, param, SymmSec, quenchtime, ll, *args):
"""
Build the Liouville operator for a system with symmetries. This
iterator returns a tuple of the row index, column index in the
Liouville space, and value of the matrix element.
**Arguments**
param : dictionary
dictionary with the setup for one simulation.
SymmSec : instance of :py:class:`exactdiag.SymmetrySector`
Allows the mapping to the basis states according to the
symmetry.
quenchtime : ``None`` or tuple of integer, float
If ``None``, access static coupling. If tuple, first entry
contains the index of the quench, the second entry contains
the time to evaluate the coupling.
ll : int
System size
pbc : bool
used to identify periodic boundary conditions. If not referenced
for the terms, \*args is used.
"""
get_op = lambda opstr : self.Operators[opstr]
dim = len(SymmSec)
# Terms of -0.5 Ld L x I and -0.5 I x (Ld L)^T
for k1, k2, elem in self.build_symm(param, SymmSec, quenchtime, ll,
*args):
for xx in range(dim):
# Ld L x I
yield (k1 * dim + xx, k2 * dim + xx, elem)
# I x (Ld L)^T (swap k1 <--> k2)
yield (dim * xx + k2, dim * xx + k1, elem)
# Terms of Li x Ld^T = Li x Li*
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Op = self.Operators[self.Op[ii]]
colidx = self._col_ind_dic(Op)
for jj in range(ll):
# ..............................................................
for k1 in range(dim):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
try:
Selk = colidx[idx1[jj]]
except KeyError:
continue
for ijk in Selk:
idx2[jj] = ijk[0]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
elem = coupl[jj] * ijk[1]
# ......................................................
for x1 in range(dim):
idx3 = SymmSec(x1)
idx4 = deepcopy(idx3)
try:
Selx = colidx[idx3[jj]]
except KeyError:
continue
for ijk2 in Selx:
idx4[jj] = ijk2[0]
try:
x2 = SymmSec[tuple(idx4)]
except KeyError:
continue
elem2 = elem * np.conj(ijk2[1])
yield (k1 * dim + x1, k2 * dim + x2, elem2)
[docs]class MBSLind(_MPOTerm):
"""
Add a many-body string Lindblad :math:`L_i` operator defined as
:math:`L_i = weight \cdot hparam \cdot \prod_{j=1}^{n} Op[j]_{i+j-1}`
acting on every sequence of :math:`n` neighboring sites. The number
of operators is determined from the length of the input list.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
the operators on the local Hilbert space used.
**Variables**
Op : list
The ii-th entry contains the n string identifiers for the
n operators of the ii-th bond term.
weight : list
The ii-th float in the list contains the weight for the
ii-th bond term.
hp : list
The ii-th string in the list contains links to the Hamiltonian
paramater with the coupling used for the ii-th Lindblad term.
**Details**
To add a MBSLind, you have to pass the following arguments to
:py:func:`mpo.MPO.AddMPOTerm`:
termtype : str
String must be equal to ``MBSLind``.
terms : list of str
String identifiers for the string of operators.
weight : float (optional keyword argument)
Overall weight for the rule
Default to 1.0.
hparam : str (optional keyword argument)
String identifying the Hamiltonian parameters allowing for quenching
the parameter or changing it for different static simulations.
Default to 1.0 (``one``).
"""
def __init__(self, Operators):
self.Operators = Operators
self.Op = []
self.weight = []
self.hp = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add a new MBSLind term. For arguments see class description
:py:class:`mpoterm.MBSLind`.
"""
self.ignored(['weight', 'hparam'], kwargs)
self.default({'weight' : 1.0, 'hparam' : 'one'}, kwargs)
if(not hasattr(terms, '__len__')):
raise Exception("terms passed into AddMPOTerm for a MBString is "
+ "not a list of Operators tags!")
if(len(terms) < 2):
raise Exception("terms passed into AddMPOTerm for a MBString "
+ "should be a list with at least two elements.")
for elem in terms:
if(not elem in self.Operators):
raise MPSOperatorNotFound((elem, "MBSLind.AddTerm"))
self.Op.append(terms)
self.weight.append(kwargs['weight'])
self.hp.append(kwargs['hparam'])
return
[docs] def write(self, fh, IntHparam, IntOperators, Params):
"""
Write the information for fortran interfaces. For description
of arguments see :py:func:`mpoterm.SiteTerm.write`.
"""
fh.write('%16i\n'%(len(self)))
for ii in range(len(self)):
hpw = 0 if(self.hp[ii] == 'one') else IntHparam[self.hp[ii]]
Hstring = '%30.15E'%(self.weight[ii])
Hstring += '%16i%16i\n'%(hpw, len(self.Op[ii]))
for jj in range(len(self.Op[ii])):
Hstring += '%16i'%(IntOperators[self.Op[ii][jj]])
fh.write(Hstring + '\n')
return
[docs] def print_term(self):
"""
Print equation for SiteLindXX terms.
"""
if(len(self) > 0):
print('MBSLind not printed in equation yet.')
return ''
[docs] def build_nosymm(self, param, sparse, quenchtime, ll, ld, *args):
"""
Build the effective hamiltonian for a system without symmetries.
This iterator returns the matrix on the complete Hilbert space
for the effective Hamiltonian. For description of terms look into
:py:func:`mpoterm.SiteTerm.build_nosymm`.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Ops = -0.5 * eye(1)
for kk in range(len(self.Op[ii])):
Opkk = get_op(self.Op[ii][kk])
Ops = kron(Ops, Opkk.conj().transpose().dot(Opkk))
jmax = ll - len(self.Op[ii]) + 1
neye = ll - len(self.Op[ii])
for jj in range(jmax):
yield kron(kron(eye(ld**jj), coupl[jj] * Ops),
eye(ld**(neye - jj)))
[docs] def lbuild_nosymm(self, param, sparse, quenchtime, ll, ld, *args):
"""
Build the Liouville operator for a system without symmetries.
This iterator returns the matrix build on the complete Liouville
space. Arguments see :py:func:`mpoterm.SiteLindXX.lbuild_nosymm`.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
dim = ld**ll
# Terms of -0.5 Ld L x I and -0.5 I x (Ld L)^T
for elem in self.build_nosymm(param, sparse, quenchtime, ll, ld):
yield kron(elem, eye(dim))
yield kron(eye(dim), elem.transpose())
# Terms of Li x Ld^T = Li x Li*
for ii in range(len(self)):
coupl = self.hparam(ii, param, quenchtime)
Ops = eye(1)
for kk in range(len(self.Op[ii])):
Opkk = get_op(self.Op[ii][kk])
Ops = kron(Ops, Opkk)
jmax = ll - len(self.Op[ii]) + 1
neye = ll - len(self.Op[ii])
for jj in range(jmax):
yield kron(kron(eye(ld**jj),
coupl[jj] * kron(Ops, eye(ld**neye))),
kron(Ops.conj(), eye(ld**(neye - jj))))
[docs] def build_symm(self, param, SymmSec, quenchtime, ll, *args):
"""
Build the effective Hamiltonian for a system with symmetries. This
iterator returns a tuple of the row index, column index, and value
of the matrix element. For complete description look into
:py:func:`mpoterm.SiteTerm.build_symm`. NOT IMPLEMENTED.
"""
for ii in range(len(self)):
# Length string
ls = len(self.Op[ii])
coupl = self.hparam(ii, param, quenchtime)
Ops = []
colidxs = []
for jj in range(ls):
Op = self.Operators[self.Op[ii][jj]]
Op = Op.conj().transpose().dot(Op)
colidx = self._col_ind_dic(Op)
Ops.append(Op)
colidxs.append(colidx)
for jj in range(ll - ls + 1):
for k1 in range(len(SymmSec)):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
Sels = []
nsel = []
try:
for kk in range(ls):
Sel = colidxs[kk][idx1[jj + kk]]
Sels.append(Sel)
nsel.append(len(Sel))
except KeyError:
continue
for selidx in self.iter_sels(nsel):
elem = -0.5 * coupl[jj]
for kk in range(ls):
idx2[jj + kk] = Sels[kk][selidx[kk]][0]
elem *= Sels[kk][selidx[kk]][1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
yield (k1, k2, elem)
[docs] def lbuild_symm(self, param, SymmSec, quenchtime, ll, *args):
"""
Build the Liouville operator for a system with symmetries. This
iterators returns a tuple of the row index, column index in the
Liouville space, and value of the matrix element. NOT IMPLEMENTED.
"""
dim = len(SymmSec)
# Terms of -0.5 Ld L x I and -0.5 I x (Ld L)^T
for k1, k2, elem in self.build_symm(param, SymmSec, quenchtime, ll,
*args):
for xx in range(dim):
# Ld L x I
yield (k1 * dim + xx, k2 * dim + xx, elem)
# I x (Ld L)^T (swap k1 <--> k2)
yield (dim * xx + k2, dim * xx + k1, elem)
# Terms of Li x Ld^T = Li * Li*
for ii in range(len(self)):
# Length string
ls = len(self.Op[ii])
coupl = self.hparam(ii, param, quenchtime)
Ops = []
colidxs = []
for jj in range(ls):
Op = self.Operators[self.Op[ii][jj]]
colidx = self._col_ind_dic(Op)
Ops.append(Op)
colidxs.append(colidx)
for jj in range(ll - ls + 1):
# ..............................................................
for k1 in range(dim):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
Sels_a = []
nsel_a = []
try:
for kk in range(ls):
Sel = colidxs[kk][idx1[jj + kk]]
Sels_a.append(Sel)
nsel_a.append(len(Sel))
except:
continue
for selidx_a in self.iter_sels(nsel_a):
elem = coupl[jj]
for kk in range(ls):
idx2[jj + kk] = Sels_a[kk][selidx_a[kk]][0]
elem *= Sels_a[kk][selidx_a[kk]][1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
# ......................................................
for x1 in range(dim):
idx3 = SymmSec(x1)
idx4 = deepcopy(idx3)
Sels_b = []
nsel_b = []
try:
for kk in range(ls):
Sel = colidxs[kk][idx3[jj + kk]]
Sels_b.append(Sel)
nsel_b.append(len(Sel))
except KeyError:
continue
for selidx_b in self.iter_sels(nsel_b):
elem2 = 1.0
for kk in range(ls):
idx4[jj + kk] = Sels_b[kk][selidx_b[kk]][0]
elem2 *= Sels_b[kk][selidx_b[kk]][1]
try:
x2 = SymmSec[tuple(idx4)]
except KeyError:
continue
elem2 = elem * np.conj(elem2)
yield (k1 * dim + x1, k2 * dim + x2, elem2)
[docs] def iter_sels(self, nsel):
"""
Iterate through selections on multiple sites.
"""
nn = np.prod(np.array(nsel))
mm = len(nsel)
selidx = [0] * mm
for ii in range(nn):
yield selidx
if(ii + 1 < nn):
selidx[0] += 1
for jj in range(1, mm):
if(selidx[jj - 1] + 1 == nsel[jj - 1]):
selidx[jj - 1] = 0
selidx[jj] += 1
else:
break
[docs]class MBSLindXY(_MPOTerm):
"""
Add a many-body string Lindblad with different arguments for L and
Ldagger.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
the operators on the local Hilbert space used.
**Variables**
Op : list
The length of the list is :math:`2 n`, where the Lindblad
operator spans n sites. The ii-th and the (n + ii)-th
string of the list define the operators acting on the
ii-th site of the Lindblad operator as L and Ldagger,
respectively.
weight : list
The ii-th float in the list contains the weight for the
ii-th bond term.
hp : list
The ii-th string in the list contains links to the Hamiltonian
paramater with the coupling used for the ii-th Lindblad term.
**Details**
To add a MBSLind, you have to pass the following arguments to
:py:func:`mpo.MPO.AddMPOTerm`:
termtype : str
String must be equal to ``MBSLindXY``.
terms : list of str
String identifiers for the string of operators.
weight : float (optional keyword argument)
Overall weight for the rule
Default to 1.0.
hparam : str (optional keyword argument)
String identifying the Hamiltonian parameters allowing for quenching
the parameter or changing it for different static simulations.
Default to 1.0 (``one``).
"""
def __init__(self, Operators):
self.Operators = Operators
self.Op = []
self.weight = []
self.hp = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add a new MBSLindXY term. For arguments see class description
:py:class:`mpoterm.MBSLindXY`.
"""
self.ignored(['weight', 'hparam'], kwargs)
self.default({'weight' : 1.0, 'hparam' : 'one'}, kwargs)
if(not hasattr(terms, '__len__')):
raise Exception("terms passed into AddMPOTerm for a MBString is "
+ "not a list of Operators tags!")
if(len(terms) < 4):
raise Exception("terms passed into AddMPOTerm for a MBString "
+ "should be a list with at least four elements.")
for elem in terms:
if(not elem in self.Operators):
raise MPSOperatorNotFound((elem, "MBSLind.AddTerm"))
self.Op.append(terms)
self.weight.append(kwargs['weight'])
self.hp.append(kwargs['hparam'])
return
[docs] def write(self, fh, IntHparam, IntOperators, Params):
"""
Write the information for fortran interfaces. For description
of arguments see :py:func:`mpoterm.SiteTerm.write`.
"""
fh.write('%16i\n'%(len(self)))
for ii in range(len(self)):
hpw = 0 if(self.hp[ii] == 'one') else IntHparam[self.hp[ii]]
Hstring = '%30.15E'%(self.weight[ii])
Hstring += '%16i%16i\n'%(hpw, len(self.Op[ii]) // 2)
for jj in range(len(self.Op[ii])):
Hstring += '%16i'%(IntOperators[self.Op[ii][jj]])
fh.write(Hstring + '\n')
return
[docs] def print_term(self):
"""
Print equation for MBSLindXY terms.
"""
if(len(self) > 0):
print('MBSLindXY not printed in equation yet.')
return ''
[docs] def build_nosymm(self, param, sparse, quenchtime, ll, ld, *args):
"""
Build the effective hamiltonian for a system without symmetries.
This iterator returns the matrix on the complete Hilbert space
for the effective Hamiltonian. For description of terms look into
:py:func:`mpoterm.SiteTerm.build_nosymm`.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
for ii in range(len(self)):
os = len(self.Op[ii]) // 2
coupl = self.hparam(ii, param, quenchtime)
Ops = -0.5 * eye(1)
for kk in range(os):
OpkkA = get_op(self.Op[ii][kk])
OpkkB = get_op(self.Op[ii][os + kk])
Ops = kron(Ops, OpkkB.conj().transpose().dot(OpkkA))
jmax = ll - os + 1
neye = ll - os
for jj in range(jmax):
yield kron(kron(eye(ld**jj), coupl[jj] * Ops),
eye(ld**(neye - jj)))
[docs] def lbuild_nosymm(self, param, sparse, quenchtime, ll, ld, *args):
"""
Build the Liouville operator for a system without symmetries.
This iterator returns the matrix build on the complete Liouville
space. Arguments see :py:func:`mpoterm.SiteLindXX.lbuild_nosymm`.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
dim = ld**ll
# Terms of -0.5 Ld L x I and -0.5 I x (Ld L)^T
for elem in self.build_nosymm(param, sparse, quenchtime, ll, ld):
yield kron(elem, eye(dim))
yield kron(eye(dim), elem.transpose())
# Terms of Li x Ld^T = Li x Li*
for ii in range(len(self)):
os = len(self.Op[ii]) // 2
coupl = self.hparam(ii, param, quenchtime)
OpsA = eye(1)
for kk in range(os):
OpkkA = get_op(self.Op[ii][kk])
OpsA = kron(Ops, OpkkA)
OpsB = eye(1)
for kk in range(os):
OpkkB = get_op(self.Op[ii][os + kk])
OpsB = kron(OpsB, OpkkB)
jmax = ll - os + 1
neye = ll - os
for jj in range(jmax):
yield kron(kron(eye(ld**jj),
coupl[jj] * kron(OpsA, eye(ld**neye))),
kron(OpsB.conj(), eye(ld**(neye - jj))))
[docs] def build_symm(self, param, SymmSec, quenchtime, ll, *args):
"""
Build the effective Hamiltonian for a system with symmetries. This
iterator returns a tuple of the row index, column index, and value
of the matrix element. For complete description look into
:py:func:`mpoterm.SiteTerm.build_symm`. NOT IMPLEMENTED.
"""
for ii in range(len(self)):
# Length string
ls = len(self.Op[ii]) // 2
coupl = self.hparam(ii, param, quenchtime)
Ops = []
colidxs = []
for jj in range(ls):
OpA = self.Operators[self.Op[ii][jj]]
OpB = self.Operators[self.Op[ii][ls + jj]]
Op = OpB.conj().transpose().dot(OpA)
colidx = self._col_ind_dic(Op)
Ops.append(Op)
colidxs.append(colidx)
for jj in range(ll - ls + 1):
for k1 in range(len(SymmSec)):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
Sels = []
nsel = []
try:
for kk in range(ls):
Sel = colidxs[kk][idx1[jj + kk]]
Sels.append(Sel)
nsel.append(len(Sel))
except KeyError:
continue
for selidx in self.iter_sels(nsel):
elem = -0.5 * coupl[jj]
for kk in range(ls):
idx2[jj + kk] = Sels[kk][selidx[kk]][0]
elem *= Sels[kk][selidx[kk]][1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
yield (k1, k2, elem)
[docs] def lbuild_symm(self, param, SymmSec, quenchtime, ll, *args):
"""
Build the Liouville operator for a system with symmetries. This
iterators returns a tuple of the row index, column index in the
Liouville space, and value of the matrix element. NOT IMPLEMENTED.
"""
dim = len(SymmSec)
# Terms of -0.5 Ld L x I and -0.5 I x (Ld L)^T
for k1, k2, elem in self.build_symm(param, SymmSec, quenchtime, ll,
*args):
for xx in range(dim):
# Ld L x I
yield (k1 * dim + xx, k2 * dim + xx, elem)
# I x (Ld L)^T (swap k1 <--> k2)
yield (dim * xx + k2, dim * xx + k1, elem)
# Terms of Li x Ld^T = Li * Li*
for ii in range(len(self)):
# Length string
ls = len(self.Op[ii]) // 2
coupl = self.hparam(ii, param, quenchtime)
colidxsA = []
colidxsB = []
for jj in range(ls):
Op = self.Operators[self.Op[ii][jj]]
colidx = self._col_ind_dic(Op)
colidxsA.append(colidx)
for jj in range(ls, 2 * ls):
Op = self.Operators[self.Op[ii][jj]]
colidx = self._col_ind_dic(Op)
colidxsB.append(colidx)
for jj in range(ll - ls + 1):
# ..............................................................
for k1 in range(dim):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
Sels_a = []
nsel_a = []
try:
for kk in range(ls):
Sel = colidxsA[kk][idx1[jj + kk]]
Sels_a.append(Sel)
nsel_a.append(len(Sel))
except:
continue
for selidx_a in self.iter_sels(nsel_a):
elem = coupl[jj]
for kk in range(ls):
idx2[jj + kk] = Sels_a[kk][selidx_a[kk]][0]
elem *= Sels_a[kk][selidx_a[kk]][1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
# ......................................................
for x1 in range(dim):
idx3 = SymmSec(x1)
idx4 = deepcopy(idx3)
Sels_b = []
nsel_b = []
try:
for kk in range(ls):
Sel = colidxsB[kk][idx3[jj + kk]]
Sels_b.append(Sel)
nsel_b.append(len(Sel))
except KeyError:
continue
for selidx_b in self.iter_sels(nsel_b):
elem2 = 1.0
for kk in range(ls):
idx4[jj + kk] = Sels_b[kk][selidx_b[kk]][0]
elem2 *= Sels_b[kk][selidx_b[kk]][1]
try:
x2 = SymmSec[tuple(idx4)]
except KeyError:
continue
elem2 = elem * np.conj(elem2)
yield (k1 * dim + x1, k2 * dim + x2, elem2)
[docs] def iter_sels(self, nsel):
"""
Iterate through selections on multiple sites.
"""
nn = np.prod(np.array(nsel))
mm = len(nsel)
selidx = [0] * mm
for ii in range(nn):
yield selidx
if(ii + 1 < nn):
selidx[0] += 1
for jj in range(1, mm):
if(selidx[jj - 1] + 1 == nsel[jj - 1]):
selidx[jj - 1] = 0
selidx[jj] += 1
else:
break
[docs]class LindExp(_MPOTerm):
"""
Add a Lindblad term with and exponential decay of the interaction in
distance between the two sites :math:`i` and :math:`j` as
:math:`|i - j|^{-1}`.
**Arguments**
**Variables**
**Details**
To add a LindExp term, you have to pass the following arguments to
:py:func:`mpo.MPO.AddMPOTerm`:
termtype : str
String must be equal to ``LindExp``
terms : list of str
String identifiers for the operators acting on the left and on the
right site.
decayparam : float (keyword argument)
Configuring the decay as 1 / decayparam**d where d is the distance.
weight : float (optional keywork argument)
Overall weight for the rule.
Default to 1.0
hparam : str (optional keyword argument)
String identifying the Hamiltonian parameters allowing for quenching
the parameter or changing it for different static simulations.
Default to 1.0 (``one``).
StringOp : str (optional keyword argument)
Possibility to add operator acting on all sites between the left and
the right site. This can be used for a Fermi phase or similar
operations (If string is equal to ``FermiPhase``, correct operators
for left and right site are constructed.)
Default to Identity (``I``).
"""
def __init__(self, Operators):
self.Operators = Operators
self.Op = []
self.weight = []
self.hp = []
self.dp = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add a new LindExp term. For arguments see class description
:py:class:`mpoterm.LindExp`.
"""
self.ignored(['weight', 'hparam', 'decayparam', 'StringOp'], kwargs)
self.required(['decayparam'], kwargs)
self.default({'weight' : 1.0, 'hparam' : 'one', 'StringOp' : 'I'},
kwargs)
weight = kwargs['weight']
hparam = kwargs['hparam']
decayp = kwargs['decayparam']
so = kwargs['StringOp']
if(not (isinstance(terms, list) or isinstance(terms, tuple))):
raise ValueError('There are no decomposition requests for'
+ 'Lindblad operators of type LindExp.')
# Local operators as product
# --------------------------
if(terms[0] not in self.Operators):
raise MPSOperatorNotFound((terms[0], "LindExp.AddTerm"))
if(terms[1] not in self.Operators):
raise MPSOperatorNotFound((terms[0], "LindExp.AddTerm"))
if(so not in self.Operators):
raise MPSOperatorNotFound((terms[0], "LindExp.AddTerm"))
if(so != 'FermiPhase'):
# No Fermi phase
# ..............
self.add([terms[0], terms[1], so], weight, hparam, decayp)
return
# Fermi phase operator
# ....................
rightname, rightwei = self.phased_operator(terms[0], 'LindExp.AddTerm')
self.add([rightname, terms[1], so], rightwei * weight, hparam, decayp)
return
[docs] def add(self, terms, weight, hparam, decayparam):
"""
Add LindExp for internal use.
**Arguments**
terms: list of three string
Two string identifiers for the operators as list.
weight: float
Weight for the whole bond term.
hparam: str
String identifier for the Hamiltonian parameter.
decayparam : float
Decay parameter of the exponential rule.
"""
self.Op.append(terms)
self.weight.append(weight)
self.hp.append(hparam)
self.dp.append(decayparam)
[docs] def write(self, fh, IntHparam, IntOperators, Params):
"""
Write the information for fortran interfaces. For description of
arguments see :py:func:`mpoterm.SiteTerm.write`.
"""
fh.write('%16i\n'%(len(self)))
for ii in range(len(self)):
hpw = 0 if(self.hp[ii] == 'one') else IntHparam[self.hp[ii]]
dp = self.dp[ii]
Hstring = '%30.15E'%(self.weight[ii])
Hstring += '%30.15E'%(Params[dp] if(isinstance(dp, str)) else dp)
Hstring += '%16i'%(hpw)
for jj in range(len(self.Op[ii])):
Hstring += '%16i'%(IntOperators[self.Op[ii][jj]])
fh.write(Hstring + '\n')
[docs] def print_term(self):
"""
Print equation for exponential terms.
"""
if(len(self) == 0): return ''
print('Warning: LindExp are not yet printed.')
return ''
[docs]class LindInfFunc(_MPOTerm):
"""
Add Lindblad operators acting on two sites where the coupling is decaying
with increasing distance between the sites.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
the operators on the local Hilbert space used.
ExpTermMPO : instance of :py:class:`mpoterm.ExpTerm`
Add the fitted exponential terms to this instance of a class.
PostProcess : bool
If ``True`` the functions are not fitted. If ``False`` fitting
procedure expresses InfiniteFunc as ExpTerms.
**Variables**
etmpo : instance of :py:class:`mpoterm.ExpTerm`
Add the fitted exponential terms to this instance of a class. Passed
to ``__init__`` as ``ExpTermMPO``.
Op : list
The ii-th entry contains the two string identifiers for the
left and right operators of the ii-th bond term.
weight : list
The ii-th float in the list contains the weight for the
ii-th bond term.
hp : list
The ii-th string in the list contains links to the Hamiltonian
paramater with the coupling used for the ii-th bond term.
yvals : list of vectors
The ii-th vector specifies the coupling with distance.
**Details**
To add a InfiniteFunc term, you have to pass the following arguments to
:py:func:`mpo.MPO.AddMPOTerm`:
termtype : str
String must be equal to ``InfiniteFunction``.
terms : list of str
strings are the name of the operators acting on all neighboring sites.
weight : float (optional keyword argument)
Overall weight for the rule
Default to 1.0.
hparam : str (optional keyword argument)
String identifying the Hamiltonian parameters allowing for quenching
the parameter or changing it for different static simulations.
Default to 1.0 (``one``).
StringOp : str (optional keyword argument)
Possibility to add operator acting on all sites between the left and
the right site. This can be used for a Fermi phase or similar
operations (If string is equal to ``FermiPhase``, correct operators
for left and right site are constructed.)
Default to Identity (``I``).
maxnterms : int (optional keyword argument)
Number of exponentials maximally allowed to fit InfiniteFunction
Default to 16.
tol : float (optional keyword argument)
This sets the tolerance for the fit of the infinite function
with the exponentials.
Default to 1e-6
func : callable (optional keyword argument)
Function returning the coupling as a function for the distance.
Combined with the argument ``L``, it may be used to specify the
InfiniteFunction
Default to ``None``
L : int (optional keyword argument)
range for fitting the exponentials to the function ``func``
Default to ``None``
x : 1d vector (optional keyword argument)
The values of the function to be fitted can be given as two vector
x and y = f(x). These are the x-values.
Default to ``None``
y : 1d vector (optional keyword argument)
The values of the function to be fitted can be given as two vector
x and y = f(x). These are values of the function.
Default to ``None``
ED : bool (optional keyword argument)
It can only be used to indicate for the EDLib not to fit an
InfiniteFunction with exponentials.
Default to ``False``
Although func, L, x, y are all marked as optional, either x and y or
func and L have to be provided.
"""
def __init__(self, Operators, ExpTermMPO, PostProcess):
self.Operators = Operators
self.etmpo = ExpTermMPO
self.PostProcess = PostProcess
self.Op = []
self.weight = []
self.hp = []
self.yvals = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add a new Lindblad infinite function terms. For arguments
see class description :py:class:`mpoterm.LindInfFunc`.
"""
# Print warnings / check required arguments / set default values
args = ['weight', 'hparam', 'StringOp', 'maxnterms', 'tol', 'func',
'L', 'x', 'y', 'ED']
defaults = {'weight' : 1.0, 'hparam' : 'one', 'Phase' : None,
'maxnterms' : 16, 'tol' : 1e-6, 'func' : None,
'L' : None, 'x' : None, 'y' : None, 'ED' : False}
self.ignored(args, kwargs)
self.default(defaults, kwargs)
if(self.PostProcess):
# Need to add dummy term for checking hparam in AddQuench
self.etmpo.AddTerm(terms, hparam=kwargs['hparam'], decayparam=1.0,
weight=0.0, Phase=kwargs['Phase'])
return
# Check x, y, func, and L directly
xx = kwargs['x']
yy = kwargs['y']
ll = kwargs['L']
func = kwargs['func']
if((xx is None) and (yy is None) and (func is not None)
and (ll is not None)):
# Generate values from function
xx = np.linspace(1, ll, ll)
yy = func(xx)
elif((xx is not None) and (yy is not None)):
pass
else:
raise Exception("Either func and L or x and y should be supplied to"
+ " AddInfiniteFunctionMPOterm!")
if(kwargs['ED']):
raise NotImplementedError('AddTermED for LindInfFunc.')
#self.AddTermED(terms, kwargs['weight'], kwargs['hparam'],
# kwargs['Phase'], xx, yy)
#return
# MPS InfiniteFunc are fitted by series of exponentials
pp = fittoexpSum(xx, yy, kwargs['maxnterms'], kwargs['tol'])
nn = len(pp) // 2
print('Number of exponentials used to fit infinite function', nn)
fitparams = []
for ii in range(nn):
fitparams.append([pp[2 * ii], pp[2 * ii + 1]])
# Local operators or decomposition request will be treated inside
# of AddexponentialMPOterm
kwargs2 = {'hparam' : kwargs['hparam'], 'StringOp' : kwargs['StringOp']}
for kk in fitparams:
kwargs2['weight'] = kwargs['weight'] * kk[0]
kwargs2['decayparam'] = kk[1]
self.etmpo.AddTerm(terms, **kwargs2)
return
[docs] def print_term(self):
"""
Print equation for InfiniteFunction terms.
"""
if(len(self) == 0): return ''
print('LindInfFunc terms not printed in equation yet.')
return ''
[docs]class LindLSpec(_MPOTerm):
"""
Add Lindblad term based on the spectral decomposition of a small
subsystem. In contrast to the :py:class:`mpoterm.LindSpec` this
approach is not based on a operator.
"""
def __init__(self, Operators):
self.Operators = Operators
self.Op = []
self.weight = []
self.hp = []
self.Te = []
self.kb = []
self.mu = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add Lindblad operators based on a spectral decomposition of
a subpart of the system. For arguments see class description of
:py:class:`mpoterm.LindLSpec`.
"""
arglist = ['weight', 'hparam', 'Te', 'kb', 'mu']
self.ignored(arglist, kwargs)
defaultlist = {'weight' : 1.0, 'hparam' : 'one', 'Te' : 0.0,
'kb' : 1.0, 'mu' : 0.0}
self.default(defaultlist, kwargs)
if(len(self) > 0):
raise Exception("Cannot add two LindLSpec at present.")
self.Op.append(terms)
self.weight.append(kwargs['weight'])
self.hp.append(kwargs['hparam'])
self.Te.append(kwargs['Te'])
self.kb.append(kwargs['kb'])
self.mu.append(kwargs['mu'])
return
[docs] def write(self, fh, IntHparam, IntOperators, Params):
"""
Not implemented.
"""
raise NotImplementedError
[docs] def print_term(self):
"""
Print equation for LindLSpec term.
"""
if(len(self) != 0):
print('LindLSpec not printed in equation yet.')
return ''
[docs] def build_nosymm(self, param, sparse, quenchtime, ll, ld, Ham):
"""
Build the effective Hamiltonian for a system without symmetries
This iterator returns the matrix on the complete Hilbert space.
**Arguments**
Ham : instance of :py:class:`mpo.MPO`
Represents the MPO of the full system. At present couplings
must be translational invariant.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
pbc = False
# Make copies / set value
locparam = deepcopy(param)
locham = deepcopy(Ham)
for ii in range(1):
coupl = self.hparam(ii, param, quenchtime)
locparam['L'] = self.Op[ii]
Trans = self.get_transitions(locparam, locham, quenchtime,
ll, ld, pbc)
for jj in range(ll - self.Op[ii] + 1):
for Top, we in Trans:
tmp = -0.5 * Top.conj().transpose().dot(Top)
yield kron(kron(eye(ld**jj), coupl[jj] * we * tmp),
eye(ld**(ll - self.Op[ii] - jj)))
[docs] def lbuild_nosymm(self, param, sparse, quenchtime, ll, ld, Ham):
"""
Build the Liouville operator for a system without symmetries
This iterator returns the matrix built on the complete Liouville
space. Arguments see :py:func:`mpoterm.SiteLindXX.lbuild_nosymm`.
Ham : instance of :py:class:`mpo.MPO`
Represents the MPO of the full system. At present couplings
must be translational invariant.
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
dim = ld**ll
pbc = False
print('>>> >>> lbuild LindLSpec')
# Terms of -0.5 Ld L x I and -0.5 I x (Ld L)^T
for elem in self.build_nosymm(param, sparse, quenchtime, ll, ld, Ham):
yield kron(elem, eye(dim))
yield kron(eye(dim), elem.transpose())
# Terms of Li Ld^T = Li x Li*
# ---------------------------
locparam = deepcopy(param)
locham = deepcopy(Ham)
for ii in range(1):
coupl = self.hparam(ii, param, quenchtime)
locparam['L'] = self.Op[ii]
Trans = self.get_transitions(locparam, locham, quenchtime,
ll, ld, pbc)
llp = ll - self.Op[ii]
for jj in range(llp + 1):
for Top, we in Trans:
yield kron(kron(eye(ld**jj),
kron(coupl[jj] * we * Top,
eye(ld**(llp)))),
kron(Top.conj(), eye(ld**(llp - jj))))
[docs] def get_transitions(self, locparam, locham, quenchtime, ll, ld, pbc):
"""
Construct the transition Lindblad operators.
"""
# Build Hamiltonian, spectrum (sparse=False)
maxdim = ld**self.Op[0] + 1
Mat = locham._build_hamiltonian_nosymm(locparam, False, maxdim,
quenchtime, pbc)
vals, vecs = la.eigh(Mat)
# Get values for thermal distribution
Te = locparam[self.Te[0]] if(isinstance(self.Te[0], str)) else self.Te[0]
kb = locparam[self.kb[0]] if(isinstance(self.kb[0], str)) else self.kb[0]
mu = locparam[self.mu[0]] if(isinstance(self.mu[0], str)) else self.mu[0]
# Select transitions (Hamiltonian must be translational invariant)
Trans = []
for k1 in range(vals.shape[0]):
for k2 in range(k1 + 1, vals.shape[0]):
# Values ascending: dE > 0
dE = vals[k2] - vals[k1]
# Do not allow transition between degenerate states
if(np.abs(dE) < 1e-13): continue
if(dE < 0.0): raise ValueError
# Calculate Fermi-Dirac probability
fh = fermi_dirac_stat(dE, Te, kb, mu)
fc = 1 - fh
Top = np.outer(np.conj(vecs[k2]), vecs[k1])
Trans.append((Top, fh))
Trans.append((Top.conj().transpose(), fc))
return Trans
[docs]class LindSpec(_MPOTerm):
"""
Add Lindblad terms for reservoirs based on the spectral
decomposition.
**Arguments**
weight : float
General weight.
Default to 1.0
hparam : str
That is the coupling of the operator to the field in
the reservoir given as
:math:`\sqrt{\\frac{1}{2 \epsilon_0 \\nu}} d_{12}`.
Each non-zero coupling means a coupling to the same
reservoir. The term is squared for translational invariant
coupling and multiplied as :math:`g[i] \cdot g[j]` if acting
on sites i and j.
Default to ``one``.
Te : str or float
Temperature of the reservoir. If str, the corresponding key
must be in the simulation dictionary.
Default to 0.0
Tapprox : tuple of (float, int, int)
Tolerance and minimal and maximal number of transitions.
If ints are None, maximal number of terms is assumed.
Default to (1e-15, None, None)
wtol : float
Tolerance for truncating a transition times
:math:`\langle e| Op | e' \\rangle`.
Default to 1e-15
ctol : float
Tolerance for truncating a whole term.
Default to 1e-15
kb : str or float
Boltzmann constant. If str, the corresponding key
must be in the simulation dictionary.
Default to 1.0
mu : str or float
Chemical potential. If str, the corresponding key
must be in the simulation dictionary.
Default to 0.0
al : str or float
Lattice constant or distance between sites. The sites are assumed
to be equidistant. If str, the corresponding key
must be in the simulation dictionary.
Default to 1.0
cc : str or float
Speed of light. If str, the corresponding key
must be in the simulation dictionary.
Default to 1.0
pvtable : bool or instance of Interpol
If False, no principal value is used and thus there is no effective
Hamiltonian. If True, principal value part is used, but no table is
given. If pvtable is not of type bool, it is assumed that an
instance of Interpol is passed.
Default to True.
"""
def __init__(self, Operators):
self.Operators = Operators
self.Op = []
self.weight = []
self.hp = []
self.table = []
self.Tapprox = []
self.wtol = []
self.ctol = []
self.idx = []
[docs] def AddTerm(self, terms, **kwargs):
"""
Add Lindblad operators based on a spectral decomposition. For
arguments see class description of :py:class:`mpoterm.LindSpec`.
"""
# First we have to set the model
self.default({'table' : CouplSpinInBEC('tmp.pkl', 1)}, kwargs)
model = kwargs.pop('table')
# Get the default list and add arguments
defaultlist = model.get_default_list()
defaultlist['weight'] = 1.0
defaultlist['hparam'] = 'one'
defaultlist['Tapprox'] = (1e-15, None, None)
defaultlist['wtol'] = 1e-15
defaultlist['ctol'] = 1e-15
# Get the arglist
arglist = list(defaultlist.keys())
# Run checks and set arguments
self.ignored(arglist, kwargs)
self.default(defaultlist, kwargs)
# Append local terms
self.Op.append(terms)
self.weight.append(kwargs['weight'])
self.hp.append(kwargs['hparam'])
self.table.append(model)
self.Tapprox.append(kwargs['Tapprox'])
self.wtol.append(kwargs['wtol'])
self.ctol.append(kwargs['ctol'])
idx = model.AddTerm(**kwargs)
self.idx.append(idx)
return
[docs] def write(self, fh, IntHparam, IntOperators, Params):
"""
Not implemented.
"""
raise NotImplementedError
[docs] def print_term(self):
"""
Print equation for LindSpec term.
"""
if(len(self) != 0):
print('LindSpec not printed in equation yet.')
return ''
[docs] def build_nosymm(self, param, sparse, quenchtime, ll, ld, Ham):
"""
Build the effective Hamiltonian for a system without symmetries.
This iterator returns the matrix on the complete Hilbert space.
**Arguments**
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
for we, L1, L2 in self.get_lindblads_nosymm(param, quenchtime, ll, ld,
Ham):
# coupl considered in get_lindblads_nosymm
#coupl = self.hparam(ii, param, quenchtime)
# Construct operator
L21 = L2.conj().transpose().dot(L1)
# Hamiltonian piece (-1j for consistency)
#yield -1j * IMAG(we) * L2.conj().transpose().dot(L1)
# Lindblad piece for effective Hamiltonian (no 0.5 necessary !!!)
#yield - REAL(we) * L2.conj().transpose().dot(L1)
yield -we * L21
[docs] def build_symm(self, param, SymmSec, quenchtime, ll, Ham):
"""
"""
for we, L1, L2 in self.get_lindblads_symm(param, SymmSec, quenchtime,
ll, Ham):
L21 = L2.conj().transpose().dot(L1)
yield -we * L21
[docs] def lbuild_nosymm(self, param, sparse, quenchtime, ll, ld, Ham):
"""
Build the Liouville operator for a system without symmetries.
This iterator returns the matrix build on the complete Liouville
space.
**Arguments**
Ham : 2d matrix
The Hamiltonian represented as matrix or sparse matrix. (Not
the rule set of the Hamiltonian.)
"""
kron, eye, get_op = self.get_sparse_dense_func(sparse, quenchtime)
dim = ld**ll
for we, L1, L2 in self.get_lindblads_nosymm(param, quenchtime, ll, ld,
Ham):
# Construct operator and yield Hamiltonian piece (do not call
# build_ham to avoid constructing all operators twice!)
L21 = L2.conj().transpose().dot(L1)
if(sparse):
L1 = sp.csc_matrix(L1)
L2 = sp.csc_matrix(L2)
L21 = sp.csc_matrix(L21)
yield kron(L21, -we * eye(dim))
yield kron(-np.conj(we) * eye(dim), L21.transpose())
# Lindblad piece L1 x L2dT (account for factor 2 here!)
wr = np.real(we)
yield kron(2 * wr * L1, L2.conj())
[docs] def lbuild_symm(self, param, SymmSec, quenchtime, ll, Ham):
"""
"""
kron, eye, get_op = self.get_sparse_dense_func(False, quenchtime)
# Dimension of the Hilbert space
spdim = len(SymmSec)
for we, L1, L2 in self.get_lindblads_symm(param, SymmSec, quenchtime,
ll, Ham):
L21 = L2.conj().transpose().dot(L1)
yield kron(L21, -we * eye(spdim))
yield kron(-np.conj(we) * eye(spdim), L21.transpose())
# Lindblad piece L1 x L2dT (account for factor 2 here!)
wr = np.real(we)
yield kron(2 * wr * L1, L2.conj())
[docs] def get_lindblads_nosymm(self, param, quenchtime, ll, ld, Ham):
"""
Generate all Lindblad operators from the spectral decomposition.
"""
for we, L1, L2 in self._get_lindblads_both(param, quenchtime, ll,
Ham, ld=ld):
yield we, L1, L2
[docs] def get_lindblads_symm(self, param, SymmSec, quenchtime, ll, Ham):
"""
Generate all Lindblad operators from the spectral decomposition.
"""
for we, L1, L2 in self._get_lindblads_both(param, quenchtime, ll,
Ham, SymmSec=SymmSec):
yield we, L1, L2
[docs] def _get_lindblads_both(self, param, quenchtime, ll, Ham, ld=None, SymmSec=None):
"""
Generate all Lindblad operators from the spectral decomposition.
"""
# Get the eigenstates to be included (search for max temperature)
ii_ = 0
Tmax = 0.0
for ii in range(len(self)):
model = self.table[ii]
Te = model[(self.idx[ii], 'Te', param)]
if(Te > Tmax):
Tmax = Te
ii_ = ii
vals, vecs, nn = self.select_eigenstates(param, Ham, ii_)
# Initialize transitions
trans = None
# Get value from parameters dictionary or as float number
get_value = lambda pa, ke : pa[ke] if(isinstance(ke, str)) else ke
for ii in range(len(self)):
# Get values of parameters
ctol = get_value(param, self.ctol[ii])
gx = self.hparam(ii, param, quenchtime)
# Model cannot be passed in parameters, we need to know that
# when adding things.
model = self.table[ii]
Te = model[(self.idx[ii], 'Te', param)]
if(not isinstance(model, (CouplSpinInBEC, CouplSpinInEM3D))):
raise Exception('Unknown coupling scheme for LindSpec.')
if(trans is None):
# Get the transitions to be included
if(SymmSec is None):
# System with no symmetry: ld must be given
trans = self.select_transitions_nosymm(param, ld, ll,
vals, vecs, nn, ii)
else:
trans = self.select_transitions_symm(param, SymmSec, ll,
vals, vecs, nn, ii)
# Precalculate results
# --------------------
_gridr = []
_omega = []
for key1 in trans.keys():
i1, j1 = key1
nk1 = len(trans[key1])
dE1 = trans[key1][0][2]
if(np.abs(dE1) < 1e-16):
# No transition for small energy differences due
# to E^3 p(E) = 0 (should be taken care of in
# select_transitions
continue
# Must skip 1 / 0
if((Te != 0.0) and np.isinf(bose_einstein_stat(np.abs(dE1), Te,
1.0, 0.0))):
continue
for key2 in trans.keys():
i2, j2 = key2
nk2 = len(trans[key2])
dE2 = trans[key2][0][2]
# Take care of delta function
if(np.abs(dE2 - dE1) > 1e-16): continue
for k1_ in range(nk1):
k1 = trans[key1][k1_][0]
w1 = trans[key1][k1_][1]
for k2_ in range(nk2):
try:
k2 = trans[key2][k2_][0]
except IndexError:
raise Exception('IndexError'
+ str(trans[key2][k2_]))
if(gx[k1] * gx[k2] == 0.0): continue
w2 = trans[key2][k2_][1]
# Calculate distance
di = np.abs(k2 - k1)
_gridr.append(di)
_omega.append(dE1)
if(len(_omega) == model.threads):
# Start calculating as soon as we reach the
# number of threads
model.pre(self.idx[ii], _gridr, _omega, param)
_gridr = []
_omega = []
if(len(_omega) != 0):
# Calculate remaining data
model.pre(self.idx[ii], _gridr, _omega, param)
# Yield the results to build the Liouville operator
# -------------------------------------------------
for key1 in trans.keys():
i1, j1 = key1
nk1 = len(trans[key1])
dE1 = trans[key1][0][2]
if(np.abs(dE1) < 1e-16):
# No transition for small energy differences due
# to E^3 p(E) = 0 (should be taken care of in
# select_transitions
continue
# Must skip 1 / 0
if((Te != 0.0) and np.isinf(bose_einstein_stat(np.abs(dE1), Te,
1.0, 0.0))):
continue
for key2 in trans.keys():
i2, j2 = key2
nk2 = len(trans[key2])
dE2 = trans[key2][0][2]
# Take care of delta function
if(np.abs(dE2 - dE1) > 1e-16): continue
for k1_ in range(nk1):
k1 = trans[key1][k1_][0]
w1 = trans[key1][k1_][1]
for k2_ in range(nk2):
try:
k2 = trans[key2][k2_][0]
except IndexError:
raise Exception('IndexError'
+ str(trans[key2][k2_]))
if(gx[k1] * gx[k2] == 0.0): continue
w2 = trans[key2][k2_][1]
# Calculate distance
di = np.abs(k2 - k1)
we = model(self.idx[ii], di, dE1, param)
# Scale with overlaps from eigenstates
we *= w1 * w2 * gx[k1] * gx[k2]
# Do not add small contributions
if(np.abs(we) < ctol): continue
# Construct Lindblads
L1 = np.outer(np.conj(vecs[:, i1]), vecs[:, j1])
L2 = np.outer(np.conj(vecs[:, i2]), vecs[:, j2])
yield we, L1, L2
if(ii + 1 < len(self)):
if(self.Op[ii] == self.Op[ii + 1]):
# Same operator, same transitions
pass
else:
# Different operator, rebuild transitions
trans = None
[docs] def select_eigenstates(self, param, Ham, ii_):
"""
Diagonalize the Hamiltonian and truncate high energy states
if desired.
"""
# Get eigenvalues
if(not isinstance(Ham, np.ndarray)):
# Sparse matrix
Ham = Ham.toarray()
vals, vecs = la.eigh(Ham)
# Generate probability distribution
model = self.table[ii_]
Te = model[(self.idx[ii_], 'Te', param)]
kb = model[(self.idx[ii_], 'kb', param)]
if(Te == 0.0):
Pd = np.zeros(vals.shape)
Pd[0] = 1.0
else:
beta = 1.0 / Te / kb
Pd = np.exp(- beta * vals)
Pd /= np.sum(Pd)
if(isinstance(self.Tapprox[0], str)):
Tapprox = param[self.Tapprox[0]]
else:
Tapprox = self.Tapprox[0]
# Get cut
nn = vals.shape[0]
for ii in range(vals.shape[0]):
if(Pd[ii] < Tapprox[0]):
nn = ii
break
if(Tapprox[1] is not None):
# Set minimum number of eigenstates
nn = max(nn, Tapprox[1])
if(Tapprox[2] is not None):
# Set maximum number of eigenstates
nn = min(nn, Tapprox[2])
# Select
vals = vals[:nn]
vecs = vecs[:, :nn]
return vals, vecs, nn
[docs] def select_transitions_nosymm(self, param, ld, ll, vals, vecs, nn, oo):
"""
Get the transitions energies between all selected eigenstates. The
returned dictionary has the a tuple of two indices of the two
corresponding eigenstates as key. The value of dictionary are two
list containing the site indices and the weight.
**Arguments**
oo : int
Choose the oo-th term of Lindblad operators.
"""
# Dimension of the Hilbert space
dim = ld**ll
# Retrieve operator
Op = self.Operators[self.Op[oo]]
# Lists for indices, sites, weight, and energy difference
transidx = []
transkk = []
transwe = []
transde = []
for ii in range(nn):
for jj in range(nn):
ket = np.reshape(vecs[:, jj], [ld] * ll)
dE = vals[ii] - vals[jj]
# Neglect zero energy transitions canceled in integral
# due to transition energy cubed.
if(np.abs(dE) < 1e-16): continue
transidx.append((ii, jj))
transde.append(dE)
transkk.append([])
transwe.append([])
for kk in range(ll):
perm = list(range(1, kk + 1)) + [0] \
+ list(range(kk + 1, ll))
tmp = np.tensordot(Op, ket, axes=([1], [kk]))
tmp = np.reshape(np.transpose(tmp, perm), [dim])
we = np.conj(vecs[:, ii]).dot(tmp)
#if(np.abs(we) > 1e-14): print('Weight for transition', we)
if(np.abs(we) >= self.wtol[0]):
transkk[-1].append((kk))
transwe[-1].append(we)
# Sort them according to energy
idx = np.argsort(np.array(transde))
# Put them into a dictionary with eigenstate indices as key
# and site, weight
transdic = {}
for ii in range(len(transidx)):
nj = len(transkk[idx[ii]])
if(nj == 0): continue
# At least one entry
transdic[transidx[idx[ii]]] = []
for jj in range(nj):
transdic[transidx[idx[ii]]].append((transkk[idx[ii]][jj],
transwe[idx[ii]][jj],
transde[idx[ii]]))
return transdic
[docs] def select_transitions_symm(self, param, SymmSec, ll, vals, vecs, nn, oo):
"""
Get the transition energies between all selected eigenstates. The
returned dicitionary has a tuple of two indices of the two
corresponding eigenstates as key. The value of the dictionary are
two lists containing the site indices and the weight.
oo : int
Choose the oo-th term of Lindblad operators.
"""
# Dimension of the Hilbert space
spdim = len(SymmSec)
# Retrieve operator
Op = self.Operators[self.Op[oo]]
colidx = self._col_ind_dic(Op)
# Lists for indices, sites, weight, and energy differences
transidx = []
transkk = []
transwe = []
transde = []
for ii in range(nn):
for jj in range(nn):
dE = vals[ii] - vals[jj]
# Neglect zero energy transitions canceled in integral
# due to transition energy cubed.
if(np.abs(dE) < 1e-16): continue
transidx.append((ii, jj))
transde.append(dE)
transkk.append([])
transwe.append([])
for kk in range(ll):
Phi = np.zeros((spdim), dtype=vecs.dtype)
# ---> Begin multiplication Phi = Op x vecs[:, jj]
for k1 in range(len(SymmSec)):
idx1 = SymmSec(k1)
idx2 = deepcopy(idx1)
skip = False
try:
Sel = colidx[idx1[kk]]
except KeyError:
skip = True
if(skip): continue
for ijk in Sel:
idx2[kk] = ijk[0]
#elem = coupl[kk] * ijk[1]
elem = 1.0 * ijk[1]
try:
k2 = SymmSec[tuple(idx2)]
except KeyError:
continue
Phi[k1] += elem * vecs[k2, jj]
# End multiplication Phi = Op x vecs[:, jj] <---
we = np.conj(vecs[:, ii]).dot(Phi)
if(np.abs(we) >= self.wtol[0]):
transkk[-1].append((kk))
transwe[-1].append(we)
# Sort them according to energy
idx = np.argsort(np.array(transde))
# Put the into a dictionary with the eigenstate indices as key
# and site, weight
transdic = {}
for ii in range(len(transidx)):
nj = len(transkk[idx[ii]])
if(nj == 0): continue
# At least one entry
transdic[transidx[idx[ii]]] = []
for jj in range(nj):
transdic[transidx[idx[ii]]].append((transkk[idx[ii]][jj],
transwe[idx[ii]][jj],
transde[idx[ii]]))
return transdic
[docs]def _append2nparray(flnm, val):
"""
Append a value to a numpy array on the hard disk. Intend for
debugging.
**Arguments**
flnm : str
Filename is existing numpy array and
will be updated with the additional entry.
val : float
Value to be appended to the numpy array.
"""
arr = np.load(flnm)
arr = list(arr)
arr.append(val)
np.save(flnm, np.array(arr))
return
[docs]class MPSOperatorNotFound(Exception):
"""
Excpetion raised if term is not found in OperatorList / dictionary
**Arguments**
ret_val : int
return code of the fortran library.
"""
def __init__(self, ret_val):
self.val = ret_val
return
[docs] def __str__(self):
"""
String representing error message.
"""
return "Operator " + str(self.val[0]) + " passed to " \
+ self.val[1] + " is not in Operator list!"