# This file is part of QuTiP: Quantum Toolbox in Python.
#
# Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# 1. Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names
# of its contributors may be used to endorse or promote products derived
# from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
###############################################################################
__all__ = ['liouvillian', 'liouvillian_ref', 'lindblad_dissipator',
'operator_to_vector', 'vector_to_operator', 'mat2vec', 'vec2mat',
'vec2mat_index', 'mat2vec_index', 'spost', 'spre', 'sprepost']
import scipy.sparse as sp
import numpy as np
from qutip.qobj import Qobj
from qutip.fastsparse import fast_csr_matrix, fast_identity
from qutip.sparse import sp_reshape
from qutip.cy.spmath import zcsr_kron
from functools import partial
[docs]def liouvillian(H, c_ops=[], data_only=False, chi=None):
"""Assembles the Liouvillian superoperator from a Hamiltonian
and a ``list`` of collapse operators. Like liouvillian, but with an
experimental implementation which avoids creating extra Qobj instances,
which can be advantageous for large systems.
Parameters
----------
H : Qobj or QobjEvo
System Hamiltonian.
c_ops : array_like of Qobj or QobjEvo
A ``list`` or ``array`` of collapse operators.
Returns
-------
L : Qobj or QobjEvo
Liouvillian superoperator.
"""
if isinstance(c_ops, (Qobj, QobjEvo)):
c_ops = [c_ops]
if chi and len(chi) != len(c_ops):
raise ValueError('chi must be a list with same length as c_ops')
h = None
if H is not None:
if isinstance(H, QobjEvo):
h = H.cte
else:
h = H
if h.isoper:
op_dims = h.dims
op_shape = h.shape
elif h.issuper:
op_dims = h.dims[0]
op_shape = [np.prod(op_dims[0]), np.prod(op_dims[0])]
else:
raise TypeError("Invalid type for Hamiltonian.")
else:
# no hamiltonian given, pick system size from a collapse operator
if isinstance(c_ops, list) and len(c_ops) > 0:
if isinstance(c_ops[0], QobjEvo):
c = c_ops[0].cte
else:
c = c_ops[0]
if c.isoper:
op_dims = c.dims
op_shape = c.shape
elif c.issuper:
op_dims = c.dims[0]
op_shape = [np.prod(op_dims[0]), np.prod(op_dims[0])]
else:
raise TypeError("Invalid type for collapse operator.")
else:
raise TypeError("Either H or c_ops must be given.")
sop_dims = [[op_dims[0], op_dims[0]], [op_dims[1], op_dims[1]]]
sop_shape = [np.prod(op_dims), np.prod(op_dims)]
spI = fast_identity(op_shape[0])
td = False
L = None
if isinstance(H, QobjEvo):
td = True
def H2L(H):
if H.isoper:
return -1.0j * (spre(H) - spost(H))
else:
return H
L = H.apply(H2L)
data = L.cte.data
elif isinstance(H, Qobj):
if H.isoper:
Ht = H.data.T
data = -1j * zcsr_kron(spI, H.data)
data += 1j * zcsr_kron(Ht, spI)
else:
data = H.data
else:
data = fast_csr_matrix(shape=(sop_shape[0], sop_shape[1]))
td_c_ops = []
for idx, c_op in enumerate(c_ops):
if isinstance(c_op, QobjEvo):
td = True
if c_op.const:
c_ = c_op.cte
elif chi:
td_c_ops.append(lindblad_dissipator(c_op, chi=chi[idx]))
continue
else:
td_c_ops.append(lindblad_dissipator(c_op))
continue
else:
c_ = c_op
if c_.issuper:
data = data + c_.data
else:
cd = c_.data.H
c = c_.data
if chi:
data = data + np.exp(1j * chi[idx]) * \
zcsr_kron(c.conj(), c)
else:
data = data + zcsr_kron(c.conj(), c)
cdc = cd * c
cdct = cdc.T
data = data - 0.5 * zcsr_kron(spI, cdc)
data = data - 0.5 * zcsr_kron(cdct, spI)
if not td:
if data_only:
return data
else:
L = Qobj()
L.dims = sop_dims
L.data = data
L.superrep = 'super'
return L
else:
if not L:
l = Qobj()
l.dims = sop_dims
l.data = data
l.superrep = 'super'
L = QobjEvo(l)
else:
L.cte.data = data
for c_op in td_c_ops:
L += c_op
return L
def liouvillian_ref(H, c_ops=[]):
"""Assembles the Liouvillian superoperator from a Hamiltonian
and a ``list`` of collapse operators.
Parameters
----------
H : qobj
System Hamiltonian.
c_ops : array_like
A ``list`` or ``array`` of collapse operators.
Returns
-------
L : qobj
Liouvillian superoperator.
"""
L = -1.0j * (spre(H) - spost(H)) if H else 0
for c in c_ops:
if c.issuper:
L += c
else:
cdc = c.dag() * c
L += spre(c) * spost(c.dag()) - 0.5 * spre(cdc) - 0.5 * spost(cdc)
return L
[docs]def lindblad_dissipator(a, b=None, data_only=False, chi=None):
"""
Lindblad dissipator (generalized) for a single pair of collapse operators
(a, b), or for a single collapse operator (a) when b is not specified:
.. math::
\\mathcal{D}[a,b]\\rho = a \\rho b^\\dagger -
\\frac{1}{2}a^\\dagger b\\rho - \\frac{1}{2}\\rho a^\\dagger b
Parameters
----------
a : Qobj or QobjEvo
Left part of collapse operator.
b : Qobj or QobjEvo (optional)
Right part of collapse operator. If not specified, b defaults to a.
Returns
-------
D : qobj, QobjEvo
Lindblad dissipator superoperator.
"""
if b is None:
b = a
ad_b = a.dag() * b
if chi:
D = spre(a) * spost(b.dag()) * np.exp(1j * chi) \
- 0.5 * spre(ad_b) - 0.5 * spost(ad_b)
else:
D = spre(a) * spost(b.dag()) - 0.5 * spre(ad_b) - 0.5 * spost(ad_b)
if isinstance(a, QobjEvo) or isinstance(b, QobjEvo):
return D
else:
return D.data if data_only else D
[docs]def operator_to_vector(op):
"""
Create a vector representation of a quantum operator given
the matrix representation.
"""
if isinstance(op, QobjEvo):
return op.apply(operator_to_vector)
q = Qobj()
q.dims = [op.dims, [1]]
q.data = sp_reshape(op.data.T, (np.prod(op.shape), 1))
return q
[docs]def vector_to_operator(op):
"""
Create a matrix representation given a quantum operator in
vector form.
"""
if isinstance(op, QobjEvo):
return op.apply(vector_to_operator)
q = Qobj()
# e.g. op.dims = [ [[rows], [cols]], [1]]
q.dims = op.dims[0]
shape = (np.prod(q.dims[0]), np.prod(q.dims[1]))
q.data = sp_reshape(op.data.T, shape[::-1]).T
return q
def mat2vec(mat):
"""
Private function reshaping matrix to vector.
"""
return mat.T.reshape(np.prod(np.shape(mat)), 1)
def vec2mat(vec, shape=None):
"""
Private function reshaping vector to matrix.
"""
if shape is None:
n = int(np.sqrt(len(vec)))
shape = (n, n)
return vec.reshape(shape[::-1]).T
def vec2mat_index(N, I):
"""
Convert a vector index to a matrix index pair that is compatible with the
vector to matrix rearrangement done by the vec2mat function.
"""
j = int(I / N)
i = I - N * j
return i, j
def mat2vec_index(N, i, j):
"""
Convert a matrix index pair to a vector index that is compatible with the
matrix to vector rearrangement done by the mat2vec function.
"""
return i + N * j
[docs]def spost(A):
"""Superoperator formed from post-multiplication by operator A
Parameters
----------
A : Qobj or QobjEvo
Quantum operator for post multiplication.
Returns
-------
super : Qobj or QobjEvo
Superoperator formed from input qauntum object.
"""
if isinstance(A, QobjEvo):
return A.apply(spost)
if not isinstance(A, Qobj):
raise TypeError('Input is not a quantum object')
if not A.isoper:
raise TypeError('Input is not a quantum operator')
S = Qobj(isherm=A.isherm, superrep='super')
S.dims = [[A.dims[0], A.dims[1]], [A.dims[0], A.dims[1]]]
S.data = zcsr_kron(A.data.T,
fast_identity(np.prod(A.shape[0])))
return S
[docs]def spre(A):
"""Superoperator formed from pre-multiplication by operator A.
Parameters
----------
A : Qobj or QobjEvo
Quantum operator for pre-multiplication.
Returns
--------
super :Qobj or QobjEvo
Superoperator formed from input quantum object.
"""
if isinstance(A, QobjEvo):
return A.apply(spre)
if not isinstance(A, Qobj):
raise TypeError('Input is not a quantum object')
if not A.isoper:
raise TypeError('Input is not a quantum operator')
S = Qobj(isherm=A.isherm, superrep='super')
S.dims = [[A.dims[0], A.dims[1]], [A.dims[0], A.dims[1]]]
S.data = zcsr_kron(fast_identity(np.prod(A.shape[1])), A.data)
return S
def _drop_projected_dims(dims):
"""
Eliminate subsystems that has been collapsed to only one state due to
a projection.
"""
return [d for d in dims if d != 1]
[docs]def sprepost(A, B):
"""Superoperator formed from pre-multiplication by operator A and post-
multiplication of operator B.
Parameters
----------
A : Qobj or QobjEvo
Quantum operator for pre-multiplication.
B : Qobj or QobjEvo
Quantum operator for post-multiplication.
Returns
--------
super : Qobj or QobjEvo
Superoperator formed from input quantum objects.
"""
if isinstance(A, QobjEvo) or isinstance(B, QobjEvo):
return spre(A) * spost(B)
else:
dims = [[_drop_projected_dims(A.dims[0]),
_drop_projected_dims(B.dims[1])],
[_drop_projected_dims(A.dims[1]),
_drop_projected_dims(B.dims[0])]]
data = zcsr_kron(B.data.T, A.data)
return Qobj(data, dims=dims, superrep='super')
from qutip.qobjevo import QobjEvo