本文整理汇总了Python中qutip.steadystate.steadystate函数的典型用法代码示例。如果您正苦于以下问题:Python steadystate函数的具体用法?Python steadystate怎么用?Python steadystate使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了steadystate函数的20个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的Python代码示例。
示例1: _correlation_me_2op_2t
def _correlation_me_2op_2t(H, rho0, tlist, taulist, c_ops, a_op, b_op,
reverse=False, args=None, options=Odeoptions()):
"""
Internal function for calculating correlation functions using the master
equation solver. See :func:`correlation` for usage.
"""
if debug:
print(inspect.stack()[0][3])
if rho0 is None:
rho0 = steadystate(H, c_ops)
elif rho0 and isket(rho0):
rho0 = ket2dm(rho0)
C_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)
rho_t_list = mesolve(
H, rho0, tlist, c_ops, [], args=args, options=options).states
if reverse:
# <A(t)B(t+tau)>
for t_idx, rho_t in enumerate(rho_t_list):
C_mat[t_idx, :] = mesolve(H, rho_t * a_op, taulist,
c_ops, [b_op], args=args,
options=options).expect[0]
else:
# <A(t+tau)B(t)>
for t_idx, rho_t in enumerate(rho_t_list):
C_mat[t_idx, :] = mesolve(H, b_op * rho_t, taulist,
c_ops, [a_op], args=args,
options=options).expect[0]
return C_mat
开发者ID:dougmcnally,项目名称:qutip,代码行数:34,代码来源:correlation.py
示例2: _correlation_me_4op_2t
def _correlation_me_4op_2t(H, rho0, tlist, taulist, c_ops,
a_op, b_op, c_op, d_op, reverse=False,
args=None, options=Odeoptions()):
"""
Calculate the four-operator two-time correlation function on the form
<A(t)B(t+tau)C(t+tau)D(t)>.
See, Gardiner, Quantum Noise, Section 5.2.1
"""
if debug:
print(inspect.stack()[0][3])
if rho0 is None:
rho0 = steadystate(H, c_ops)
elif rho0 and isket(rho0):
rho0 = ket2dm(rho0)
C_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)
rho_t = mesolve(
H, rho0, tlist, c_ops, [], args=args, options=options).states
for t_idx, rho in enumerate(rho_t):
C_mat[t_idx, :] = mesolve(H, d_op * rho * a_op, taulist,
c_ops, [b_op * c_op],
args=args, options=options).expect[0]
return C_mat
开发者ID:dougmcnally,项目名称:qutip,代码行数:29,代码来源:correlation.py
示例3: _correlation_es_2op_1t
def _correlation_es_2op_1t(H, rho0, tlist, c_ops, a_op, b_op, reverse=False,
args=None, options=Odeoptions()):
"""
Internal function for calculating correlation functions using the
exponential series solver. See :func:`correlation_ss` usage.
"""
if debug:
print(inspect.stack()[0][3])
# contruct the Liouvillian
L = liouvillian(H, c_ops)
# find the steady state
if rho0 is None:
rho0 = steadystate(L)
elif rho0 and isket(rho0):
rho0 = ket2dm(rho0)
# evaluate the correlation function
if reverse:
# <A(t)B(t+tau)>
solC_tau = ode2es(L, rho0 * a_op)
return esval(expect(b_op, solC_tau), tlist)
else:
# default: <A(t+tau)B(t)>
solC_tau = ode2es(L, b_op * rho0)
return esval(expect(a_op, solC_tau), tlist)
开发者ID:dougmcnally,项目名称:qutip,代码行数:28,代码来源:correlation.py
示例4: _spectrum_es
def _spectrum_es(H, wlist, c_ops, a_op, b_op):
"""
Internal function for calculating the spectrum of the correlation function
:math:`\left<A(\\tau)B(0)\\right>`.
"""
if debug:
print(inspect.stack()[0][3])
# construct the Liouvillian
L = liouvillian(H, c_ops)
# find the steady state density matrix and a_op and b_op expecation values
rho0 = steadystate(L)
a_op_ss = expect(a_op, rho0)
b_op_ss = expect(b_op, rho0)
# eseries solution for (b * rho0)(t)
es = ode2es(L, b_op * rho0)
# correlation
corr_es = expect(a_op, es)
# covariance
cov_es = corr_es - np.real(np.conjugate(a_op_ss) * b_op_ss)
# spectrum
spectrum = esspec(cov_es, wlist)
return spectrum
开发者ID:bcriger,项目名称:qutip,代码行数:30,代码来源:correlation.py
示例5: _correlation_es_2t
def _correlation_es_2t(H, state0, tlist, taulist, c_ops, a_op, b_op, c_op):
"""
Internal function for calculating the three-operator two-time
correlation function:
<A(t)B(t+tau)C(t)>
using an exponential series solver.
"""
# the solvers only work for positive time differences and the correlators
# require positive tau
if state0 is None:
rho0 = steadystate(H, c_ops)
tlist = [0]
elif isket(state0):
rho0 = ket2dm(state0)
else:
rho0 = state0
if debug:
print(inspect.stack()[0][3])
# contruct the Liouvillian
L = liouvillian(H, c_ops)
corr_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)
solES_t = ode2es(L, rho0)
# evaluate the correlation function
for t_idx in range(len(tlist)):
rho_t = esval(solES_t, [tlist[t_idx]])
solES_tau = ode2es(L, c_op * rho_t * a_op)
corr_mat[t_idx, :] = esval(expect(b_op, solES_tau), taulist)
return corr_mat
开发者ID:JonathanUlm,项目名称:qutip,代码行数:34,代码来源:correlation.py
示例6: _spectrum_es
def _spectrum_es(H, wlist, c_ops, a_op, b_op):
"""
Internal function for calculating the spectrum of the correlation function
:math:`\left<A(\\tau)B(0)\\right>`.
"""
if debug:
print(inspect.stack()[0][3])
# construct the Liouvillian
L = liouvillian(H, c_ops)
# find the steady state density matrix and a_op and b_op expecation values
rho0 = steadystate(L)
a_op_ss = expect(a_op, rho0)
b_op_ss = expect(b_op, rho0)
# eseries solution for (b * rho0)(t)
es = ode2es(L, b_op * rho0)
# correlation
corr_es = expect(a_op, es)
# covariance
cov_es = corr_es - a_op_ss * b_op_ss
# tidy up covariance (to combine, e.g., zero-frequency components that cancel)
cov_es.tidyup()
# spectrum
spectrum = esspec(cov_es, wlist)
return spectrum
开发者ID:JonathanUlm,项目名称:qutip,代码行数:32,代码来源:correlation.py
示例7: coherence_function_g1
def coherence_function_g1(H, taulist, c_ops, a_op, solver="me", args=None,
options=Options(ntraj=[20, 100])):
"""
Calculate the normalized first-order quantum coherence function:
.. math::
g^{(1)}(\\tau) = \lim_{t \to \infty}
\\frac{\\langle a^\\dagger(t+\\tau)a(t)\\rangle}
{\\langle a^\\dagger(t)a(t)\\rangle}
using the quantum regression theorem and the evolution solver indicated by
the `solver` parameter. Note: g1 is only defined for stationary
statistics (uses steady state).
Parameters
----------
H : :class:`qutip.qobj.Qobj`
system Hamiltonian.
taulist : *list* / *array*
list of times for :math:`\\tau`. taulist must be positive and contain
the element `0`.
c_ops : list of :class:`qutip.qobj.Qobj`
list of collapse operators.
a_op : :class:`qutip.qobj.Qobj`
The annihilation operator of the mode.
solver : str
choice of solver (`me` for master-equation and
`es` for exponential series)
options : :class:`qutip.solver.Options`
solver options class. `ntraj` is taken as a two-element list because
the `mc` correlator calls `mcsolve()` recursively; by default,
`ntraj=[20, 100]`. `mc_corr_eps` prevents divide-by-zero errors in
the `mc` correlator; by default, `mc_corr_eps=1e-10`.
Returns
-------
g1: *array*
The normalized first-order coherence function.
"""
# first calculate the steady state photon number
rho0 = steadystate(H, c_ops)
n = np.array([expect(rho0, a_op.dag() * a_op)])
# calculate the correlation function G1 and normalize with n to obtain g1
G1 = correlation_2op_1t(H, None, taulist, c_ops, a_op.dag(), a_op,
args=args, solver=solver, options=options)
g1 = G1 / n
return g1
开发者ID:JonathanUlm,项目名称:qutip,代码行数:59,代码来源:correlation.py
示例8: floquet_master_equation_steadystate
def floquet_master_equation_steadystate(H, A):
"""
Returns the steadystate density matrix (in the floquet basis!) for the
Floquet-Markov master equation.
"""
c_ops = floquet_collapse_operators(A)
rho_ss = steadystate(H, c_ops)
return rho_ss
开发者ID:Marata459,项目名称:qutip,代码行数:8,代码来源:floquet.py
示例9: coherence_function_g2
def coherence_function_g2(H, rho0, taulist, c_ops, a_op, solver="me",
args=None, options=Odeoptions()):
"""
Calculate the second-order quantum coherence function:
.. math::
g^{(2)}(\\tau) =
\\frac{\\langle a^\\dagger(0)a^\\dagger(\\tau)a(\\tau)a(0)\\rangle}
{\\langle a^\\dagger(\\tau)a(\\tau)\\rangle
\\langle a^\\dagger(0)a(0)\\rangle}
Parameters
----------
H : :class:`qutip.qobj.Qobj`
system Hamiltonian.
rho0 : :class:`qutip.qobj.Qobj`
Initial state density matrix (or state vector). If 'rho0' is
'None', then the steady state will be used as initial state.
taulist : *list* / *array*
list of times for :math:`\\tau`.
c_ops : list of :class:`qutip.qobj.Qobj`
list of collapse operators.
a_op : :class:`qutip.qobj.Qobj`
The annihilation operator of the mode.
solver : str
choice of solver (currently only 'me')
Returns
-------
g2, G2: tuble of *array*
The normalized and unnormalized second-order coherence function.
"""
# first calculate the photon number
if rho0 is None:
rho0 = steadystate(H, c_ops)
n = np.array([expect(rho0, a_op.dag() * a_op)])
else:
n = mesolve(
H, rho0, taulist, c_ops, [a_op.dag() * a_op],
args=args, options=options).expect[0]
# calculate the correlation function G2 and normalize with n to obtain g2
G2 = correlation_4op_1t(H, rho0, taulist, c_ops,
a_op.dag(), a_op.dag(), a_op, a_op,
solver=solver, args=args, options=options)
g2 = G2 / (n[0] * n)
return g2, G2
开发者ID:dougmcnally,项目名称:qutip,代码行数:58,代码来源:correlation.py
示例10: countstat_current
def countstat_current(L, c_ops=None, rhoss=None, J_ops=None):
"""
Calculate the current corresponding a system Liouvillian `L` and a list of
current collapse operators `c_ops` or current superoperators `J_ops`
(either must be specified). Optionally the steadystate density matrix
`rhoss` and a list of current superoperators `J_ops` can be specified. If
either of these are omitted they are computed internally.
Parameters
----------
L : :class:`qutip.Qobj`
Qobj representing the system Liouvillian.
c_ops : array / list (optional)
List of current collapse operators.
rhoss : :class:`qutip.Qobj` (optional)
The steadystate density matrix corresponding the system Liouvillian
`L`.
J_ops : array / list (optional)
List of current superoperators.
Returns
--------
I : array
The currents `I` corresponding to each current collapse operator
`c_ops` (or, equivalently, each current superopeator `J_ops`).
"""
if J_ops is None:
if c_ops is None:
raise ValueError("c_ops must be given if J_ops is not")
J_ops = [sprepost(c, c.dag()) for c in c_ops]
if rhoss is None:
if c_ops is None:
raise ValueError("c_ops must be given if rhoss is not")
rhoss = steadystate(L, c_ops)
rhoss_vec = mat2vec(rhoss.full()).ravel()
N = len(J_ops)
I = np.zeros(N)
for i, Ji in enumerate(J_ops):
I[i] = expect_rho_vec(Ji.data, rhoss_vec, 1)
return I
开发者ID:NunoEdgarGub1,项目名称:qutip,代码行数:50,代码来源:countstat.py
示例11: _correlation_me_2t
def _correlation_me_2t(H, state0, tlist, taulist, c_ops, a_op, b_op, c_op,
args={}, options=Options()):
"""
Internal function for calculating the three-operator two-time
correlation function:
<A(t)B(t+tau)C(t)>
using a master equation solver.
"""
# the solvers only work for positive time differences and the correlators
# require positive tau
if state0 is None:
rho0 = steadystate(H, c_ops)
tlist = [0]
elif isket(state0):
rho0 = ket2dm(state0)
else:
rho0 = state0
if debug:
print(inspect.stack()[0][3])
rho_t = mesolve(H, rho0, tlist, c_ops, [],
args=args, options=options).states
corr_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)
H_shifted, c_ops_shifted, _args = _transform_L_t_shift(H, c_ops, args)
if config.tdname:
_cython_build_cleanup(config.tdname)
rhs_clear()
for t_idx, rho in enumerate(rho_t):
if not isinstance(H, Qobj):
_args["_t0"] = tlist[t_idx]
corr_mat[t_idx, :] = mesolve(
H_shifted, c_op * rho * a_op, taulist, c_ops_shifted,
[b_op], args=_args, options=options
).expect[0]
if t_idx == 1:
options.rhs_reuse = True
if config.tdname:
_cython_build_cleanup(config.tdname)
rhs_clear()
return corr_mat
开发者ID:ajgpitch,项目名称:qutip,代码行数:47,代码来源:correlation.py
示例12: _correlation_me_4op_1t
def _correlation_me_4op_1t(H, rho0, tlist, c_ops, a_op, b_op, c_op, d_op, args=None, options=Options()):
"""
Calculate the four-operator two-time correlation function on the form
<A(0)B(tau)C(tau)D(0)>.
See, Gardiner, Quantum Noise, Section 5.2.1
"""
if debug:
print(inspect.stack()[0][3])
if rho0 is None:
rho0 = steadystate(H, c_ops)
elif rho0 and isket(rho0):
rho0 = ket2dm(rho0)
return mesolve(H, d_op * rho0 * a_op, tlist, c_ops, [b_op * c_op], args=args, options=options).expect[0]
开发者ID:ntezak,项目名称:qutip,代码行数:17,代码来源:correlation.py
示例13: _correlation_me_2op_1t
def _correlation_me_2op_1t(H, rho0, tlist, c_ops, a_op, b_op, reverse=False, args=None, options=Options()):
"""
Internal function for calculating correlation functions using the master
equation solver. See :func:`correlation_ss` for usage.
"""
if debug:
print(inspect.stack()[0][3])
if rho0 is None:
rho0 = steadystate(H, c_ops)
elif rho0 and isket(rho0):
rho0 = ket2dm(rho0)
if reverse:
# <A(t)B(t+tau)>
return mesolve(H, rho0 * a_op, tlist, c_ops, [b_op], args=args, options=options).expect[0]
else:
# <A(t+tau)B(t)>
return mesolve(H, b_op * rho0, tlist, c_ops, [a_op], args=args, options=options).expect[0]
开发者ID:ntezak,项目名称:qutip,代码行数:20,代码来源:correlation.py
示例14: _correlation_es_2op_2t
def _correlation_es_2op_2t(H, rho0, tlist, taulist, c_ops, a_op, b_op,
reverse=False, args=None, options=Odeoptions()):
"""
Internal function for calculating correlation functions using the
exponential series solver. See :func:`correlation` usage.
"""
if debug:
print(inspect.stack()[0][3])
# contruct the Liouvillian
L = liouvillian(H, c_ops)
if rho0 is None:
rho0 = steadystate(L)
elif rho0 and isket(rho0):
rho0 = ket2dm(rho0)
C_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)
solES_t = ode2es(L, rho0)
# evaluate the correlation function
if reverse:
# <A(t)B(t+tau)>
for t_idx in range(len(tlist)):
rho_t = esval(solES_t, [tlist[t_idx]])
solES_tau = ode2es(L, rho_t * a_op)
C_mat[t_idx, :] = esval(expect(b_op, solES_tau), taulist)
else:
# default: <A(t+tau)B(t)>
for t_idx in range(len(tlist)):
rho_t = esval(solES_t, [tlist[t_idx]])
solES_tau = ode2es(L, b_op * rho_t)
C_mat[t_idx, :] = esval(expect(a_op, solES_tau), taulist)
return C_mat
开发者ID:dougmcnally,项目名称:qutip,代码行数:38,代码来源:correlation.py
示例15: _spectrum_pi
def _spectrum_pi(H, wlist, c_ops, a_op, b_op, use_pinv=False):
"""
Internal function for calculating the spectrum of the correlation function
:math:`\left<A(\\tau)B(0)\\right>`.
"""
L = H if issuper(H) else liouvillian(H, c_ops)
tr_mat = tensor([qeye(n) for n in L.dims[0][0]])
N = np.prod(L.dims[0][0])
A = L.full()
b = spre(b_op).full()
a = spre(a_op).full()
tr_vec = np.transpose(mat2vec(tr_mat.full()))
rho_ss = steadystate(L)
rho = np.transpose(mat2vec(rho_ss.full()))
I = np.identity(N * N)
P = np.kron(np.transpose(rho), tr_vec)
Q = I - P
spectrum = np.zeros(len(wlist))
for idx, w in enumerate(wlist):
if use_pinv:
MMR = np.linalg.pinv(-1.0j * w * I + A)
else:
MMR = np.dot(Q, np.linalg.solve(-1.0j * w * I + A, Q))
s = np.dot(tr_vec,
np.dot(a, np.dot(MMR, np.dot(b, np.transpose(rho)))))
spectrum[idx] = -2 * np.real(s[0, 0])
return spectrum
开发者ID:JonathanUlm,项目名称:qutip,代码行数:37,代码来源:correlation.py
示例16: countstat_current_noise
def countstat_current_noise(L, c_ops, rhoss=None, J_ops=None, R=False):
"""
Compute the cross-current noise spectrum for a list of collapse operators
`c_ops` corresponding to monitored currents, given the system
Liouvillian `L`. The current collapse operators `c_ops` should be part
of the dissipative processes in `L`, but the `c_ops` given here does not
necessarily need to be all collapse operators contributing to dissipation
in the Liouvillian. Optionally, the steadystate density matrix `rhoss`
and/or the pseudo inverse `R` of the Liouvillian `L`, and the current
operators `J_ops` correpsonding to the current collapse operators `c_ops`
can also be specified. If `R` is not given, the cross-current correlations
will be computed directly without computing `R` explicitly. If either of
`rhoss` and `J_ops` are omitted, they will be computed internally.
Parameters
----------
L : :class:`qutip.Qobj`
Qobj representing the system Liouvillian.
c_ops : array / list
List of current collapse operators.
rhoss : :class:`qutip.Qobj` (optional)
The steadystate density matrix corresponding the system Liouvillian
`L`.
J_ops : array / list (optional)
List of current superoperators.
R : :class:`qutip.Qobj` (optional)
Qobj representing the pseudo inverse of the system Liouvillian `L`.
Returns
--------
I, S : tuple of arrays
The currents `I` corresponding to each current collapse operator
`c_ops` (or, equivalently, each current superopeator `J_ops`) and the
zero-frequency cross-current correlation `S`.
"""
if rhoss is None:
rhoss = steadystate(L, c_ops)
if J_ops is None:
J_ops = [sprepost(c, c.dag()) for c in c_ops]
rhoss_vec = mat2vec(rhoss.full()).ravel()
N = len(J_ops)
I = np.zeros(N)
S = np.zeros((N, N))
if R:
if R is True:
R = pseudo_inverse(L, rhoss)
for i, Ji in enumerate(J_ops):
for j, Jj in enumerate(J_ops):
if i == j:
I[i] = expect_rho_vec(Ji.data, rhoss_vec, 1)
S[i, j] = I[i]
S[i, j] -= expect_rho_vec((Ji * R * Jj + Jj * R * Ji).data,
rhoss_vec, 1)
else:
N = np.prod(L.dims[0][0])
rhoss_vec = operator_to_vector(rhoss)
tr_op = tensor([identity(n) for n in L.dims[0][0]])
tr_op_vec = operator_to_vector(tr_op)
Pop = sp.kron(rhoss_vec.data, tr_op_vec.data.T, format='csc')
Iop = sp.eye(N*N, N*N, format='csc')
Q = Iop - Pop
A = L.data.tocsc()
rhoss_vec = mat2vec(rhoss.full()).ravel()
for j, Jj in enumerate(J_ops):
Qj = Q * Jj.data * rhoss_vec
X_rho_vec = sp.linalg.splu(A, permc_spec='COLAMD').solve(Qj)
for i, Ji in enumerate(J_ops):
if i == j:
S[i, i] = I[i] = expect_rho_vec(Ji.data, rhoss_vec, 1)
S[i, j] -= expect_rho_vec(Ji.data * Q, X_rho_vec, 1)
return I, S
开发者ID:JonathanUlm,项目名称:qutip,代码行数:89,代码来源:countstat.py
示例17: coherence_function_g2
def coherence_function_g2(H, state0, taulist, c_ops, a_op, solver="me", args={},
options=Options(ntraj=[20, 100])):
"""
Calculate the normalized second-order quantum coherence function:
.. math::
g^{(2)}(\\tau) =
\\frac{\\langle A^\\dagger(0)A^\\dagger(\\tau)A(\\tau)A(0)\\rangle}
{\\langle A^\\dagger(\\tau)A(\\tau)\\rangle
\\langle A^\\dagger(0)A(0)\\rangle}
using the quantum regression theorem and the evolution solver indicated by
the `solver` parameter.
Parameters
----------
H : Qobj
system Hamiltonian, may be time-dependent for solver choice of `me` or
`mc`.
state0 : Qobj
Initial state density matrix :math:`\\rho(t_0)` or state vector
:math:`\\psi(t_0)`. If 'state0' is 'None', then the steady state will
be used as the initial state. The 'steady-state' is only implemented
for the `me` and `es` solvers.
taulist : array_like
list of times for :math:`\\tau`. taulist must be positive and contain
the element `0`.
c_ops : list
list of collapse operators, may be time-dependent for solver choice of
`me` or `mc`.
a_op : Qobj
operator A.
solver : str
choice of solver (`me` for master-equation and
`es` for exponential series).
options : Options
solver options class. `ntraj` is taken as a two-element list because
the `mc` correlator calls `mcsolve()` recursively; by default,
`ntraj=[20, 100]`. `mc_corr_eps` prevents divide-by-zero errors in
the `mc` correlator; by default, `mc_corr_eps=1e-10`.
Returns
-------
g2, G2 : tuple
The normalized and unnormalized second-order coherence function.
"""
# first calculate the photon number
if state0 is None:
state0 = steadystate(H, c_ops)
n = np.array([expect(state0, a_op.dag() * a_op)])
else:
n = mesolve(H, state0, taulist, c_ops, [a_op.dag() * a_op]).expect[0]
# calculate the correlation function G2 and normalize with n to obtain g2
G2 = correlation_3op_1t(H, state0, taulist, c_ops,
a_op.dag(), a_op.dag()*a_op, a_op,
solver=solver, args=args, options=options)
g2 = G2 / (n[0] * n)
return g2, G2
开发者ID:nwlambert,项目名称:qutip,代码行数:63,代码来源:correlation.py
示例18: countstat_current_noise
def countstat_current_noise(L, c_ops, wlist=None, rhoss=None, J_ops=None,
sparse=True, method='direct'):
"""
Compute the cross-current noise spectrum for a list of collapse operators
`c_ops` corresponding to monitored currents, given the system
Liouvillian `L`. The current collapse operators `c_ops` should be part
of the dissipative processes in `L`, but the `c_ops` given here does not
necessarily need to be all collapse operators contributing to dissipation
in the Liouvillian. Optionally, the steadystate density matrix `rhoss`
and the current operators `J_ops` correpsonding to the current collapse
operators `c_ops` can also be specified. If either of
`rhoss` and `J_ops` are omitted, they will be computed internally.
'wlist' is an optional list of frequencies at which to evaluate the noise
spectrum.
Note:
The default method is a direct solution using dense matrices, as sparse
matrix methods fail for some examples of small systems.
For larger systems it is reccomended to use the sparse solver
with the direct method, as it avoids explicit calculation of the
pseudo-inverse, as described in page 67 of "Electrons in nanostructures"
C. Flindt, PhD Thesis, available online:
http://orbit.dtu.dk/fedora/objects/orbit:82314/datastreams/file_4732600/content
Parameters
----------
L : :class:`qutip.Qobj`
Qobj representing the system Liouvillian.
c_ops : array / list
List of current collapse operators.
rhoss : :class:`qutip.Qobj` (optional)
The steadystate density matrix corresponding the system Liouvillian
`L`.
wlist : array / list (optional)
List of frequencies at which to evaluate (if none are given, evaluates
at zero frequency)
J_ops : array / list (optional)
List of current superoperators.
sparse : bool
Flag that indicates whether to use sparse or dense matrix methods when
computing the pseudo inverse. Default is false, as sparse solvers
can fail for small systems. For larger systems the sparse solvers
are reccomended.
Returns
--------
I, S : tuple of arrays
The currents `I` corresponding to each current collapse operator
`c_ops` (or, equivalently, each current superopeator `J_ops`) and the
zero-frequency cross-current correlation `S`.
"""
if rhoss is None:
rhoss = steadystate(L, c_ops)
if J_ops is None:
J_ops = [sprepost(c, c.dag()) for c in c_ops]
N = len(J_ops)
I = np.zeros(N)
if wlist is None:
S = np.zeros((N, N,1))
wlist=[0.]
else:
S = np.zeros((N, N,len(wlist)))
if sparse == False:
rhoss_vec = mat2vec(rhoss.full()).ravel()
for k,w in enumerate(wlist):
R = pseudo_inverse(L, rhoss=rhoss, w= w, sparse = sparse, method=method)
for i, Ji in enumerate(J_ops):
for j, Jj in enumerate(J_ops):
if i == j:
I[i] = expect_rho_vec(Ji.data, rhoss_vec, 1)
S[i, j,k] = I[i]
S[i, j,k] -= expect_rho_vec((Ji * R * Jj
+ Jj * R * Ji).data,
rhoss_vec, 1)
else:
if method == "direct":
N = np.prod(L.dims[0][0])
rhoss_vec = operator_to_vector(rhoss)
tr_op = tensor([identity(n) for n in L.dims[0][0]])
tr_op_vec = operator_to_vector(tr_op)
Pop = sp.kron(rhoss_vec.data, tr_op_vec.data.T, format='csr')
Iop = sp.eye(N*N, N*N, format='csr')
Q = Iop - Pop
#.........这里部分代码省略.........
开发者ID:NunoEdgarGub1,项目名称:qutip,代码行数:101,代码来源:countstat.py
示例19: spectrum_ss
def spectrum_ss(H, wlist, c_ops, a_op, b_op):
"""
Calculate the spectrum corresponding to a correlation function
:math:`\left<A(\\tau)B(0)\\right>`, i.e., the Fourier transform of the
correlation function:
.. math::
S(\omega) = \int_{-\infty}^{\infty} \left<A(\\tau)B(0)\\right>
e^{-i\omega\\tau} d\\tau.
Parameters
----------
H : :class:`qutip.qobj`
system Hamiltonian.
wlist : *list* / *array*
list of frequencies for :math:`\\omega`.
c_ops : list of :class:`qutip.qobj`
list of collapse operators.
a_op : :class:`qutip.qobj`
operator A.
b_op : :class:`qutip.qobj`
operator B.
Returns
-------
spectrum: *array*
An *array* with spectrum :math:`S(\omega)` for the frequencies
specified in `wlist`.
"""
if debug:
print(inspect.stack()[0][3])
# contruct the Liouvillian
L = liouvillian(H, c_ops)
# find the steady state density matrix and a_op and b_op expecation values
rho0 = steadystate(L)
a_op_ss = expect(a_op, rho0)
b_op_ss = expect(b_op, rho0)
# eseries solution for (b * rho0)(t)
es = ode2es(L, b_op * rho0)
# correlation
corr_es = expect(a_op, es)
# covarience
cov_es = corr_es - np.real(np.conjugate(a_op_ss) * b_op_ss)
# spectrum
spectrum = esspec(cov_es, wlist)
return spectrum
开发者ID:dougmcnally,项目名称:qutip,代码行数:63,代码来源:correlation.py
示例20: spectrum_pi
def spectrum_pi(H, wlist, c_ops, a_op, b_op, use_pinv=False):
"""
Calculate the spectrum corresponding to a correlation function
:math:`\left<A(\\tau)B(0)\\right>`, i.e., the Fourier transform of the
correlation function:
.. math::
S(\omega) = \int_{-\infty}^{\infty} \left<A(\\tau)B(0)\\right>
e^{-i\omega\\tau} d\\tau.
Parameters
----------
H : :class:`qutip.qobj`
system Hamiltonian.
wlist : *list* / *array*
list of frequencies for :math:`\\omega`.
c_ops : list of :class:`qutip.qobj`
list of collapse operators.
a_op : :class:`qutip.qobj`
operator A.
b_op : :class:`qutip.qobj`
operator B.
Returns
-------
s_vec: *array*
An *array* with spectrum :math:`S(\omega)` for the frequencies
specified in `wlist`.
"""
L = H if issuper(H) else liouvillian_fast(H, c_ops)
tr_mat = tensor([qeye(n) for n in L.dims[0][0]])
N = prod(L.dims[0][0])
A = L.full()
b = spre(b_op).full()
a = spre(a_op).full()
tr_vec = transpose(mat2vec(tr_mat.full()))
rho_ss = steadystate(L)
rho = transpose(mat2vec(rho_ss.full()))
I = np.identity(N * N)
P = np.kron(transpose(rho), tr_vec)
Q = I - P
s_vec = np.zeros(len(wlist))
for idx, w in enumerate(wlist):
if use_pinv:
MMR = numpy.linalg.pinv(-1.0j * w * I + A)
else:
MMR = np.dot(Q, np.linalg.solve(-1.0j * w * I + A, Q))
s = np.dot(tr_vec, np.dot(a, np.dot(MMR, np.dot(b, transpose(rho)))))
s_vec[idx] = -2 * np.real(s[0, 0])
return s_vec
开发者ID:dougmcnally,项目名称:qutip,代码行数:70,代码来源:correlation.py
注:本文中的qutip.steadystate.steadystate函数示例由纯净天空整理自Github/MSDocs等源码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。 |
请发表评论