"""
The models constructed in this module provide a the most simplistic entry
to openMPS. For designing more sophisticated Hamiltonians, it is necessary
to go through the steps of creating the operators and building an MPO on
your own.
**Details**
Pre-built models return the operators, the MPO, and an observable class. The
following pre-built models exist.
* Quantum Ising model :py:class:`model.QuantumIsingModel`
* Bose-Hubbard model :py:class:`model.BoseHubbardModel`
* Bilinear-Biquadratic spin-one model: :py:class:`BilinearBiquadraticModel`
"""
#import MPSPyLib as mps
import numpy as np
__all__ = ['QuantumIsingModel', 'BoseHubbardModel',
'BilinearBiquadraticModel']
[docs]class _Model():
"""
Abstract class for any model providing basic functions.
"""
[docs] def build_observables(self):
"""
Build an observable class for the bilinear biquadratic model.
"""
self.Obs = mps.Observables(self.Ops)
return
[docs] def build(self):
"""
Return the operators, the MPO, and a measurement object.
"""
self.build_operators()
self.build_mpo()
self.build_observables()
return self.Ops, self.Ham, self.Obs
[docs] def __call__(self):
"""
Use the call attribute to return the classes for operators, MPOs,
and measurements.
"""
return self.build()
[docs]class QuantumIsingModel(_Model):
"""
Build the Hamiltonian including the operators and observables class
for the quantum Ising model. After creating an
instance of the model, you can call the object like a function
and it will return a tuple with the operators, the MPO, and
the observables class. The Hamiltonian is defined as:
.. math::
H = -J \sum_{i=1}^{L-1} \sigma_{i}^{z} \sigma_{i+1}^{z}
-g \sum_{i=1}^{L} \sigma_{i}^{x}
**Arguments**
No arguments
**Details**
We choose the Pauli operators to be diagonal in :math:`\sigma^{x}`
which allows to use the :math:`Z_{2}` symmetry in the Ising model.
The following string identifiers are defined for each term in the
Hamiltonian:
+-------------------+-----------------------------------------------------+
| String identifier | Hamiltonian term |
+===================+=====================================================+
| ``J`` | :math:`\sigma_{i}^{z} \sigma_{i+1}^{z}` |
+-------------------+-----------------------------------------------------+
| ``g`` | :math:`\sigma_{i}^{x}` |
+-------------------+-----------------------------------------------------+
The following operators are built based on the spin operators:
+--------+-----------------------------------------------------+
| String | Definition |
+========+=====================================================+
| I | Identity :math:`1` |
+--------+-----------------------------------------------------+
| splus | Spin raising :math:`S^{+}` (in x-direction) |
+--------+-----------------------------------------------------+
| sminus | Spin lowering :math:`S^{-}` (in x-direction) |
+--------+-----------------------------------------------------+
| sz | Pauli-z matrix :math:`\sigma^{z}` (not diagonal) |
+--------+-----------------------------------------------------+
| sx | :math:`S^{x} = 0.5 (S^{+} + S^{-})` (diagonal) |
+--------+-----------------------------------------------------+
| gen | Generator for :math:`Z_2` (0 and 1 on the diagonal) |
+--------+-----------------------------------------------------+
"""
def __init__(self):
pass
[docs] def build_operators(self):
"""
Build all necessary operators for the quantum Ising Hamiltonian.
"""
Ops = mps.BuildSpinOperators(spin=0.5)
Ops['sx'] = 2 * Ops['sz']
Ops['sz'] = - (Ops['splus'] + Ops['sminus'])
Ops['gen'] = np.array([[0, 0], [0, 1.]])
self.Ops = Ops
return
[docs] def build_mpo(self):
"""
Build the MPO for the quantum Ising model.
"""
Ham = mps.MPO(self.Ops)
# External field
Ham.AddMPOTerm('site', 'sx', hparam='g', weight=-1.0)
# Interaction
Ham.AddMPOTerm('bond', ['sz', 'sz'], hparam='J', weight=-1.0)
self.Ham = Ham
return
[docs]class XXXModel(_Model):
"""
Build the Hamiltonian including the operators and observable class
for the spin-1/2 XXX model. After creating an
instance of the model, you can call the object like a function
and it will return a tuple with the operators, the MPO, and
the observables class. The Hamiltonian is defined as:
.. math::
H = -J \sum_{i=1}^{L-1} \left( \sigma_{i}^{x} \sigma_{i+1}^{x}
+ \sigma_{i}^{y} \sigma_{i+1}^{y}
+ \sigma_{i}^{z} \sigma_{i+1}^{z} \\right)
- g \sum_{i=1}^{L} \sigma_{i}^{z}
**Arguments**
No arguments.
**Details**
The following string identifiers are defined for each term in the
Hamiltonian:
+-------------------+-----------------------------------------------------+
| String identifier | Hamiltonian term |
+===================+=====================================================+
| ``J`` | Interaction, i.e., xx-term, yy-term, and zz-term |
| | at equal strength. |
+-------------------+-----------------------------------------------------+
| ``g`` | :math:`\sigma_{i}^{z}` |
+-------------------+-----------------------------------------------------+
The following operators are built based on the spin operators:
+--------+-----------------------------------------------------+
| String | Definition |
+========+=====================================================+
| I | Identity :math:`1` |
+--------+-----------------------------------------------------+
| splus | Spin raising :math:`S^{+}` (in x-direction) |
+--------+-----------------------------------------------------+
| sminus | Spin lowering :math:`S^{-}` (in x-direction) |
+--------+-----------------------------------------------------+
| sz | Spin-z matrix :math:`\sigma^{z}` |
+--------+-----------------------------------------------------+
"""
def __init__(self):
pass
[docs] def build_operators(self):
"""
Build all necessary operators for the XXX model.
"""
Ops = mps.BuildSpinOperators(spin=0.5)
self.Ops = Ops
return
[docs] def build_mpo(self):
"""
Build the MPO for the XXX model.
"""
Ham = mps.MPO(self.Ops)
# External field
Ham.AddMPOTerm('site', 'sz', hparam='g', weight=-1.0)
# Interaction
Ham.AddMPOTerm('bond', ['splus', 'sminus'], hparam='J', weight=-0.5)
Ham.AddMPOTerm('bond', ['sz', 'sz'], hparam='J', weight=-1.0)
self.Ham = Ham
return
[docs]class XXZModel(_Model):
"""
Build the Hamiltonian including the operators and observable class
for the spin-1/2 XXZ model. After creating an
instance of the model, you can call the object like a function
and it will return a tuple with the operators, the MPO, and
the observables class. The Hamiltonian is defined as:
.. math::
H = -J \sum_{i=1}^{L-1} \left( \sigma_{i}^{x} \sigma_{i+1}^{x}
+ \sigma_{i}^{y} \sigma_{i+1}^{y} \\right)
- J_{z} \sum_{i=1}^{L-1} \sigma_{i}^{z} \sigma_{i+1}^{z}
- g \sum_{i=1}^{L} \sigma_{i}^{z}
**Arguments**
No arguments.
**Details**
The following string identifiers are defined for each term in the
Hamiltonian:
+-------------------+-----------------------------------------------------+
| String identifier | Hamiltonian term |
+===================+=====================================================+
| ``J`` | Interaction of xx-term and yy-term |
| | at equal strength. |
+-------------------+-----------------------------------------------------+
| ``Jz`` | Interaction of zz-term. |
+-------------------+-----------------------------------------------------+
| ``g`` | :math:`\sigma_{i}^{z}` |
+-------------------+-----------------------------------------------------+
The following operators are built based on the spin operators:
+--------+-----------------------------------------------------+
| String | Definition |
+========+=====================================================+
| I | Identity :math:`1` |
+--------+-----------------------------------------------------+
| splus | Spin raising :math:`S^{+}` (in x-direction) |
+--------+-----------------------------------------------------+
| sminus | Spin lowering :math:`S^{-}` (in x-direction) |
+--------+-----------------------------------------------------+
| sz | Spin-z matrix :math:`\sigma^{z}` |
+--------+-----------------------------------------------------+
"""
def __init__(self):
pass
[docs] def build_operators(self):
"""
Build all necessary operators for the XXZ model.
"""
Ops = mps.BuildSpinOperators(spin=0.5)
self.Ops = Ops
return
[docs] def build_mpo(self):
"""
Build the MPO for the XXZ model.
"""
Ham = mps.MPO(self.Ops)
# External field
Ham.AddMPOTerm('site', 'sz', hparam='g', weight=-1.0)
# Interaction
Ham.AddMPOTerm('bond', ['splus', 'sminus'], hparam='J', weight=-0.5)
Ham.AddMPOTerm('bond', ['sz', 'sz'], hparam='Jz', weight=-1.0)
self.Ham = Ham
return
[docs]class BoseHubbardModel(_Model):
"""
Build the Hamiltonian including the operators and observables class
for the Bose-Hubbard model. After creating an
instance of the model, you can call the object like a function
and it will return a tuple with the operators, the MPO, and
the observables class. The Hamiltonian is defined as:
.. math::
H = - J \sum_{i=1}^{L-1} (b_{i}^{\dagger} b_{i+1} + h.c.)
+ \\frac{U}{2} \sum_{i=1}^{L} (n_{i}^2 - n_{i})
+ \mu \sum_{i=1}^{L} n_{i}
**Arguments**
nmax : int, optional
Maximal number of bosons on one site. (If you need spinful or
flavorful bosons, minimal number of bosons, etc., you have to
go beyond the default models.
Default to 4.
chem_pot : boolean, optional
Specifies if the chemical potential term :math:`\mu` should
be included in the Hamiltonian.
Default to ``False``.
**Details**
The string identifiers for each Hamiltonian parameters are:
+-------------------+-----------------------------------------------------+
| String identifier | Hamiltonian term |
+===================+=====================================================+
| ``J`` | math:`b_{i}^{\dagger} b_{i+1} + h.c.` (Tunneling) |
+-------------------+-----------------------------------------------------+
| ``U`` | math:`n_{i}^2 - n_{i}` (On-site interaction) |
+-------------------+-----------------------------------------------------+
| ``mu`` | math:`n_{i}` (chemical potential) |
+-------------------+-----------------------------------------------------+
The following operators are built:
+-------------+----------------------------------------------+
| String | Definition |
+=============+==============================================+
| I | Identity :math:`1` |
+-------------+----------------------------------------------+
| bdagger | Creation operator :math:`b^{\dagger}` |
+-------------+----------------------------------------------+
| b | Annihilation operator :math:`b` |
+-------------+----------------------------------------------+
| nbtotal | Number operator :math:`n = b^{\dagger} b` |
+-------------+----------------------------------------------+
| interaction | :math:`0.5 (n^2 - n)` |
+-------------+----------------------------------------------+
"""
def __init__(self, nmax=4, chem_pot=False):
self.chem_pot = chem_pot
self.nmax = nmax
[docs] def build_operators(self):
"""
Build all necessary operators of the Bose-Hubbard Hamiltonian.
"""
Ops = mps.BuildBoseOperators(self.nmax)
Ops['interaction'] = 0.5 * (np.dot(Ops['nbtotal'], Ops['nbtotal'])
- Ops['nbtotal'])
self.Ops = Ops
[docs] def build_mpo(self):
"""
Build all necessary operators for the Bose-Hubbard model.
"""
Ham = mps.MPO(self.Ops)
# Tunneling
Ham.AddMPOTerm('bond', ['bdagger','b'], hparam='J', weight=-1.0)
# On-site interaction
Ham.AddMPOTerm('site', 'interaction', hparam='U')
if(self.chem_pot):
Ham.AddMPOTerm('site', 'nbtotal', param='mu')
self.Ham = Ham
[docs]class BilinearBiquadraticModel(_Model):
"""
Build the Hamiltonian including the operators and observables class
for the bilinear biquadratic model. After creating an
instance of the model, you can call the object like a function
and it will return a tuple with the operators, the MPO, and
the observables class. The Hamiltonian is defined as:
.. math::
H = J \sum_{i=1}^{L-1} \left( \cos(\\theta) \\vec{S}_{i} \\vec{S}_{i+1}
+ \sin(\\theta) \left( \\vec{S}_{i} \\vec{S}_{i+1} \\right)^2 \\right)
- D \sum_{i=1}^{L} \left( S_{i}^{z} \\right)^2
- \delta S_{1}^{z}
**Arguments**
symmetry_breaking : bool, optional
Specifies if a symmetry breaking field should be included
in the Hamiltonian.
Default to ``False``.
fix_theta : float or None, optional
If ``None`` (default), the more general Hamiltonian with
thirteen bond rules is built. If you know that theta is
fix for all simulations (statics and dynamics), you may
pass the actual value reducing the problem to nine bond
rules.
use_symm : bool, optional
We can use the symmetry of the Bilinear biquadratic model. This
requires that theta is given. A warning is printed that symmetry
is not used if theta is not given.
**Details**
The string identifiers for each Hamiltonian parameter are:
+-------------------+------------------------------------------------------+
| String identifier | Hamiltonian term |
+===================+======================================================+
| ``cost`` | math:`\\vec{S}_{i} \\vec{S}_{i+1}` |
+-------------------+------------------------------------------------------+
| ``sint`` | math:`\left( \\vec{S}_{i} \\vec{S}_{i+1} \\right)^2` |
+-------------------+------------------------------------------------------+
| ``D`` | math:`\left( S_{i}^{z} \\right)^2` |
+-------------------+------------------------------------------------------+
| ``Delta`` | math:`S_{1}^{z}` |
+-------------------+------------------------------------------------------+
The following operators are built.
+--------+---------------------------------------+-------------------------+
| String | Definition | Comments |
+========+=======================================+=========================+
| I | Identity :math:`1` | - |
+--------+---------------------------------------+-------------------------+
| splus | Spin raising :math:`S^{+}` | - |
+--------+---------------------------------------+-------------------------+
| sminus | Spin lowering :math:`S^{-}` | - |
+--------+---------------------------------------+-------------------------+
| sz | Spin z :math:`S^{z}` | - |
+--------+---------------------------------------+-------------------------+
| sx | :math:`S^{x} = 0.5 (S^{+} + S^{-})` | Only without symmetry |
+--------+---------------------------------------+-------------------------+
| sy1 | :math:`S^{y1} = 0.5 (S^{+} - S^{-})` | Only without symmetry |
+--------+---------------------------------------+-------------------------+
| sy2 | :math:`S^{y2} = 0.5 (-S^{+} + S^{-})` | Only without symmetry |
+--------+---------------------------------------+-------------------------+
| sx2 | :math:`S^{x^2}` | Only without symmetry |
+--------+---------------------------------------+-------------------------+
| sy2 | :math:`S^{y^2}` | Only without symmetry |
+--------+---------------------------------------+-------------------------+
| sz2 | :math:`S^{z^2}` | - |
+--------+---------------------------------------+-------------------------+
| sxz | :math:`S^{x} S^{z}` | Only without symmetry |
+--------+---------------------------------------+-------------------------+
| sxy1 | :math:`S^{x} S^{y1}` | Only without symmetry |
+--------+---------------------------------------+-------------------------+
| sxy2 | :math:`S^{x} S^{y2}` | Only without symmetry |
+--------+---------------------------------------+-------------------------+
| szy1 | :math:`S^{z} S^{y1}` | Only without symmetry |
+--------+---------------------------------------+-------------------------+
| szy2 | :math:`S^{z} S^{y2}` | Only without symmetry |
+--------+---------------------------------------+-------------------------+
"""
def __init__(self, symmetry_breaking=False, fix_theta=None, use_symm=False):
self.symmetry_breaking = symmetry_breaking
self.fix_theta = fix_theta
self.use_symm = use_symm
if((fix_theta is None) and use_symm):
print('Warning: Symmetry will not be used due to missing value '
+ 'for theta.')
self.use_symm = False
[docs] def build_operators(self):
"""
Build all necessary operators for the spin-one bilinear
biquadratic model.
"""
Ops = mps.BuildSpinOperators(1.0)
# Build Sx, Sy
Ops['sx'] = 0.5 * (Ops['splus'] + Ops['sminus'])
Ops['sya'] = 0.5 * (Ops['splus'] - Ops['sminus'])
Ops['syb'] = Ops['sya'].transpose()
# Build squared operators
Ops['sx2'] = Ops['sx'].dot(Ops['sx'])
Ops['sy2'] = Ops['sya'].dot(Ops['sya'])
Ops['sz2'] = Ops['sz'].dot(Ops['sz'])
# Build mixed operators
Ops['sxz'] = Ops['sx'].dot(Ops['sz'])
Ops['sxya'] = Ops['sx'].dot(Ops['sya'])
Ops['sxyb'] = Ops['sx'].dot(Ops['syb'])
Ops['szya'] = Ops['sz'].dot(Ops['sya'])
Ops['szyb'] = Ops['sz'].dot(Ops['syb'])
self.Ops = Ops
return
[docs] def build_mpo(self):
"""
Build the MPO for the spin-one bilinear biquadratic model.
"""
Ham = mps.MPO(self.Ops)
# External field
Ham.AddMPOTerm('site', 'sz2', hparam='D', weight=-1.0)
# Perturbation
if(self.symmetry_breaking):
Ham.AddMPOTerm('site', 'sz', hparam='Delta', weight=-1.0)
if(self.fix_theta is None):
# Build complete rule set
# .......................
# cosine terms
Ham.AddMPOTerm('bond', ['sx', 'sx'], hparam='cost')
Ham.AddMPOTerm('bond', ['sya', 'syb'], hparam='cost', weight=0.5)
Ham.AddMPOTerm('bond', ['sz', 'sz'], hparam='cost')
# sine terms
Ham.AddMPOTerm('bond', ['sx2', 'sx2'], hparam='sint')
Ham.AddMPOTerm('bond', ['sy2', 'sy2'], hparam='sint')
Ham.AddMPOTerm('bond', ['sz2', 'sz2'], hparam='sint')
Ham.AddMPOTerm('bond', ['sxz', 'sxz'], hparam='sint')
Ham.AddMPOTerm('bond', ['sxya', 'sxyb'], hparam='sint')
Ham.AddMPOTerm('bond', ['szya', 'szyb'], hparam='sint')
else:
# Build two site operator for decomposition
# -----------------------------------------
dim = self.Ops['I'].shape[0]**2
Mat = np.zeros((dim, dim))
# cosine terms
Mat += np.cos(self.fix_theta) \
* np.kron(self.Ops['sx'], self.Ops['sx'])
Mat += 0.5 * np.cos(self.fix_theta) \
* np.kron(self.Ops['sya'], self.Ops['syb'])
Mat += 0.5 * np.cos(self.fix_theta) \
* np.kron(self.Ops['syb'], self.Ops['sya'])
Mat += np.cos(self.fix_theta) \
* np.kron(self.Ops['sz'], self.Ops['sz'])
# sine terms
Mat += np.sin(self.fix_theta) \
* np.kron(self.Ops['sx2'], self.Ops['sx2'])
Mat += np.sin(self.fix_theta) \
* np.kron(self.Ops['sy2'], self.Ops['sy2'])
Mat += np.sin(self.fix_theta) \
* np.kron(self.Ops['sz2'], self.Ops['sz2'])
Mat += np.sin(self.fix_theta) \
* np.kron(self.Ops['sxz'], self.Ops['sxz'])
Mat += np.sin(self.fix_theta) \
* np.kron(self.Ops['sxz'].transpose(),
self.Ops['sxz'].transpose())
Mat += np.sin(self.fix_theta) \
* np.kron(self.Ops['sxya'], self.Ops['sxyb'])
Mat += np.sin(self.fix_theta) \
* np.kron(self.Ops['sxya'].transpose(),
self.Ops['sxyb'].transpose())
Mat += np.sin(self.fix_theta) \
* np.kron(self.Ops['szya'], self.Ops['szyb'])
Mat += np.sin(self.fix_theta) \
* np.kron(self.Ops['szya'].transpose(),
self.Ops['szyb'].transpose())
if(self.use_symm):
gen1 = np.array([0, 1, 2])
eye1 = np.diag(np.eye(3))
gen2 = (np.kron(gen1, eye1) + np.kron(eye1, gen1))%3
# Have to delete all operators not fulfilling the symmetry
del self.Ops['sx']
del self.Ops['sya']
del self.Ops['syb']
del self.Ops['sx2']
del self.Ops['sy2']
del self.Ops['sxz']
del self.Ops['sxya']
del self.Ops['sxyb']
del self.Ops['szya']
del self.Ops['szyb']
Ham.AddMPOTerm('bond', Mat, gen2=[gen2])
raise NotImplementedError("Symmetry for BLBQ not supported "
+ "in F90 subroutines.")
else:
Ham.AddMPOTerm('bond', Mat)
self.Ham = Ham
return
[docs] def get_symmetry_breaking_field(self, value, ll):
"""
Get coupling for the symmetry breaking field.
**Arguments**
value : float
Value of the symmetry breaking field
ll : int
Number of sites
"""
tmp = np.zeros((ll))
tmp[0] = value
return tmp