本文整理汇总了Python中qutip.steadystate函数的典型用法代码示例。如果您正苦于以下问题:Python steadystate函数的具体用法?Python steadystate怎么用?Python steadystate使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了steadystate函数的19个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于我们的系统推荐出更棒的Python代码示例。
示例1: _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:thomasaref,项目名称:TA_software,代码行数:30,代码来源:qutip_spectrum_look.py
示例2: test_ho_lgmres
def test_ho_lgmres():
"Steady state: Thermal HO - iterative-lgmres solver"
# thermal steadystate of an oscillator: compare numerics with analytical
# formula
a = destroy(40)
H = 0.5 * 2 * np.pi * a.dag() * a
gamma1 = 0.05
wth_vec = np.linspace(0.1, 3, 20)
p_ss = np.zeros(np.shape(wth_vec))
for idx, wth in enumerate(wth_vec):
n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature
c_op_list = []
rate = gamma1 * (1 + n_th)
c_op_list.append(np.sqrt(rate) * a)
rate = gamma1 * n_th
c_op_list.append(np.sqrt(rate) * a.dag())
rho_ss = steadystate(H, c_op_list, method='iterative-lgmres')
p_ss[idx] = np.real(expect(a.dag() * a, rho_ss))
p_ss_analytic = 1.0 / (np.exp(1.0 / wth_vec) - 1)
delta = sum(abs(p_ss_analytic - p_ss))
assert_equal(delta < 1e-3, True)
开发者ID:jrjohansson,项目名称:qutip,代码行数:25,代码来源:test_steadystate.py
示例3: test_qubit_power
def test_qubit_power():
"Steady state: Thermal qubit - power solver"
# thermal steadystate of a qubit: compare numerics with analytical formula
sz = sigmaz()
sm = destroy(2)
H = 0.5 * 2 * np.pi * sz
gamma1 = 0.05
wth_vec = np.linspace(0.1, 3, 20)
p_ss = np.zeros(np.shape(wth_vec))
for idx, wth in enumerate(wth_vec):
n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature
c_op_list = []
rate = gamma1 * (1 + n_th)
c_op_list.append(np.sqrt(rate) * sm)
rate = gamma1 * n_th
c_op_list.append(np.sqrt(rate) * sm.dag())
rho_ss = steadystate(H, c_op_list, method='power')
p_ss[idx] = expect(sm.dag() * sm, rho_ss)
p_ss_analytic = np.exp(-1.0 / wth_vec) / (1 + np.exp(-1.0 / wth_vec))
delta = sum(abs(p_ss_analytic - p_ss))
assert_equal(delta < 1e-5, True)
开发者ID:jrjohansson,项目名称:qutip,代码行数:26,代码来源:test_steadystate.py
示例4: steady
def steady(N = 20): # number of basis states to consider
n=num(N)
a = destroy(N)
H = a.dag() * a
print H.eigenstates()
#psi0 = basis(N, 10) # initial state
kappa = 0.1 # coupling to oscillator
c_op_list = []
n_th_a = 2 # temperature with average of 2 excitations
rate = kappa * (1 + n_th_a)
c_op_list.append(sqrt(rate) * a) # decay operators
rate = kappa * n_th_a
c_op_list.append(sqrt(rate) * a.dag()) # excitation operators
final_state = steadystate(H, c_op_list)
fexpt = expect(a.dag() * a, final_state)
#tlist = linspace(0, 100, 100)
#mcdata = mcsolve(H, psi0, tlist, c_op_list, [a.dag() * a], ntraj=100)
#medata = mesolve(H, psi0, tlist, c_op_list, [a.dag() * a])
#plot(tlist, mcdata.expect[0],
#plt.plot(tlist, medata.expect[0], lw=2)
plt.axhline(y=fexpt, color='r', lw=1.5) # ss expt. value as horiz line (= 2)
plt.ylim([0, 10])
plt.show()
开发者ID:priyanka27s,项目名称:TA_software,代码行数:26,代码来源:test_qutip.py
示例5: find_expect
def find_expect(self, vg, pwr, fd):
if self.power_plot:
phi, pwr=vg
else:
phi, fd=vg
pwr_fridge=pwr-self.atten
lin_pwr=0.001*10**(pwr_fridge/10.0)
Omega=sqrt(lin_pwr/h*2.0)
gamma, Delta=self._get_GammaDelta(fd=fd, f0=self.f0, Np=self.Np, gamma=self.gamma)
g_el=self._get_Gamma_C(fd=fd)
wTvec=self._get_fTvec(phi=phi, gamma=gamma, Delta=Delta, fd=fd, Psaw=lin_pwr)
if self.acoustic_plot:
Om=Omega*sqrt(gamma/fd)
else:
Om=Omega*sqrt(g_el/fd)
wT = wTvec-fd*self.nvec #rotating frame of gate drive \omega_m-m*\omega_\gate
transmon_levels = Qobj(diag(wT[range(self.N_dim)]))
rate1 = (gamma+g_el)*(1.0 + self.N_gamma)
rate2 = (gamma+g_el)*self.N_gamma
c_ops=[sqrt(rate1)*self.a_op, sqrt(rate2)*self.a_dag]#, sqrt(rate3)*self.a_op, sqrt(rate4)*self.a_dag]
Omega_vec=-0.5j*(Om*self.a_dag - conj(Om)*self.a_op)
H=transmon_levels +Omega_vec
final_state = steadystate(H, c_ops) #solve master equation
fexpt=expect(self.a_op, final_state) #expectation value of relaxation operator
#return fexpt
if self.acoustic_plot:
return 1.0*gamma/Om*fexpt
else:
return 1.0*sqrt(g_el*gamma)/Om*fexpt
开发者ID:thomasaref,项目名称:TA_software,代码行数:31,代码来源:qubit_saturation_Lamb.py
示例6: rhos
def rhos(self, nslice=None):
"""rho
return steadystate density matrices"""
self.precalc = True
if self.noisy:
print("N = {}, len() = {}".format(self.N_field_levels, len(self.long_range)))
def progress(*args):
if self.noisy:
print("|", sep="", end="", flush=True)
return args
if nslice is not None:
return np.asarray(
[progress(qt.steadystate(ham, self._c_ops())) for ham in list(self.hamiltonian())[nslice]]
).T
else:
return np.asarray([progress(qt.steadystate(ham, self._c_ops())) for ham in self.hamiltonian()]).T
开发者ID:fergusbarratt,项目名称:masters-project,代码行数:19,代码来源:quantumoptics.py
示例7: test_dqd_current
def test_dqd_current():
"Counting statistics: current and current noise in a DQD model"
G = 0
L = 1
R = 2
sz = projection(3, L, L) - projection(3, R, R)
sx = projection(3, L, R) + projection(3, R, L)
sR = projection(3, G, R)
sL = projection(3, G, L)
w0 = 1
tc = 0.6 * w0
GammaR = 0.0075 * w0
GammaL = 0.0075 * w0
nth = 0.00
eps_vec = np.linspace(-1.5*w0, 1.5*w0, 20)
J_ops = [GammaR * sprepost(sR, sR.dag())]
c_ops = [np.sqrt(GammaR * (1 + nth)) * sR,
np.sqrt(GammaR * (nth)) * sR.dag(),
np.sqrt(GammaL * (nth)) * sL,
np.sqrt(GammaL * (1 + nth)) * sL.dag(),
]
I = np.zeros(len(eps_vec))
S = np.zeros(len(eps_vec))
for n, eps in enumerate(eps_vec):
H = (eps/2 * sz + tc * sx)
L = liouvillian(H, c_ops)
rhoss = steadystate(L)
I[n], S[n] = countstat_current_noise(L, [], rhoss=rhoss, J_ops=J_ops)
I2 = countstat_current(L, rhoss=rhoss, J_ops=J_ops)
assert_(abs(I[n] - I2) < 1e-8)
I2 = countstat_current(L, c_ops, J_ops=J_ops)
assert_(abs(I[n] - I2) < 1e-8)
Iref = tc**2 * GammaR / (tc**2 * (2 + GammaR/GammaL) +
GammaR**2/4 + eps_vec**2)
Sref = 1 * Iref * (
1 - 8 * GammaL * tc**2 *
(4 * eps_vec**2 * (GammaR - GammaL) +
GammaR * (3 * GammaL * GammaR + GammaR**2 + 8*tc**2)) /
(4 * tc**2 * (2 * GammaL + GammaR) + GammaL * GammaR**2 +
4 * eps_vec**2 * GammaL)**2
)
assert_allclose(I, Iref, 1e-4)
assert_allclose(S, Sref, 1e-4)
开发者ID:NunoEdgarGub1,项目名称:qutip,代码行数:54,代码来源:test_countstat.py
示例8: find_expect
def find_expect(vg): #phi=0.1, Omega_vec=3.0):
phi, Omega=vg#.shape
Omega_vec=- 0.5j*(Omega*adag - conj(Omega)*a)
Ej = Ejmax*absolute(cos(pi*phi)) #Josephson energy as function of Phi.
wTvec = -Ej + sqrt(8.0*Ej*Ec)*(nvec+0.5)+Ecvec #\omega_m
wT = wTvec-wdvec #rotating frame of gate drive \omega_m-m*\omega_\gate
transmon_levels = Qobj(diag(wT[range(N)]))
H=transmon_levels +Omega_vec #- 0.5j*(Omega_true*adag - conj(Omega_true)*a)
final_state = steadystate(H, c_op_list) #solve master equation
return expect( a, final_state) #expectation value of relaxation operator
开发者ID:thomasaref,项目名称:TA_software,代码行数:14,代码来源:qutip_steadystate.py
示例9: find_expect
def find_expect(vg, self=a): #phi=0.1, Omega_vec=3.0):
phi, Omega=vg#.shape
Omega_vec=-0.5j*(Omega*self.a_dag - conj(Omega)*self.a_op)
Ej = self.Ejmax*absolute(cos(pi*phi)) #Josephson energy as function of Phi.
wTvec = (-Ej + sqrt(8.0*Ej*self.Ec)*(self.nvec+0.5)+self.Ecvec)/h #\omega_m
wT = wTvec-a.fdvec #rotating frame of gate drive \omega_m-m*\omega_\gate
transmon_levels = Qobj(diag(wT[range(self.N_dim)]))
H=transmon_levels +Omega_vec #- 0.5j*(Omega_true*adag - conj(Omega_true)*a)
final_state = steadystate(H, self.c_ops) #solve master equation
return expect( self.a_op, final_state) #expectation value of relaxation operator
开发者ID:thomasaref,项目名称:TA_software,代码行数:14,代码来源:qubit_saturation.py
示例10: test_driven_cavity_power_bicgstab
def test_driven_cavity_power_bicgstab():
"Steady state: Driven cavity - power-bicgstab solver"
N = 30
Omega = 0.01 * 2 * np.pi
Gamma = 0.05
a = destroy(N)
H = Omega * (a.dag() + a)
c_ops = [np.sqrt(Gamma) * a]
M = build_preconditioner(H, c_ops, method='power')
rho_ss = steadystate(H, c_ops, method='power-bicgstab', M=M, use_precond=1)
rho_ss_analytic = coherent_dm(N, -1.0j * (Omega)/(Gamma/2))
assert_((rho_ss - rho_ss_analytic).norm() < 1e-4)
开发者ID:ajgpitch,项目名称:qutip,代码行数:14,代码来源:test_steadystate.py
示例11: test_driven_cavity_lgmres
def test_driven_cavity_lgmres():
"Steady state: Driven cavity - iterative-lgmres solver"
N = 30
Omega = 0.01 * 2 * np.pi
Gamma = 0.05
a = destroy(N)
H = Omega * (a.dag() + a)
c_ops = [np.sqrt(Gamma) * a]
rho_ss = steadystate(H, c_ops, method='iterative-lgmres')
rho_ss_analytic = coherent_dm(N, -1.0j * (Omega)/(Gamma/2))
assert_((rho_ss - rho_ss_analytic).norm() < 1e-4)
开发者ID:jrjohansson,项目名称:qutip,代码行数:15,代码来源:test_steadystate.py
示例12: jc_steadystate
def jc_steadystate(self, N, wc, wa, g, kappa, gamma,
pump, psi0, use_rwa, tlist):
# Hamiltonian
a = tensor(destroy(N), identity(2))
sm = tensor(identity(N), destroy(2))
if use_rwa:
# use the rotating wave approxiation
H = wc * a.dag(
) * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
else:
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (
a.dag() + a) * (sm + sm.dag())
# collapse operators
c_op_list = []
n_th_a = 0.0 # zero temperature
rate = kappa * (1 + n_th_a)
c_op_list.append(np.sqrt(rate) * a)
rate = kappa * n_th_a
if rate > 0.0:
c_op_list.append(np.sqrt(rate) * a.dag())
rate = gamma
if rate > 0.0:
c_op_list.append(np.sqrt(rate) * sm)
rate = pump
if rate > 0.0:
c_op_list.append(np.sqrt(rate) * sm.dag())
# find the steady state
rho_ss = steadystate(H, c_op_list)
return expect(a.dag() * a, rho_ss), expect(sm.dag() * sm, rho_ss)
开发者ID:JonathanUlm,项目名称:qutip,代码行数:39,代码来源:test_mesolve.py
示例13: _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>`.
"""
#print issuper(H)
L = H if issuper(H) else liouvillian(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 = identity(N * N)
P = kron(transpose(rho), tr_vec)
Q = I - P
spectrum = zeros(len(wlist))
for idx, w in enumerate(wlist):
if use_pinv:
MMR = pinv(-1.0j * w * I + A)
else:
MMR = dot(Q, solve(-1.0j * w * I + A, Q))
s = dot(tr_vec,
dot(a, dot(MMR, dot(b, transpose(rho)))))
spectrum[idx] = -2 * real(s[0, 0])
return spectrum
开发者ID:thomasaref,项目名称:TA_software,代码行数:39,代码来源:qutip_time_test.py
示例14: find_expect
def find_expect(phi=0.1, Omega_vec=3.0, wd=wd):
wdvec=nvec*wd
Ej = Ejmax*absolute(cos(pi*phi)) #; % Josephson energy as function of Phi.
wq0=sqrt(8*Ej*Ec)
X=Np*pi*(wq0-f0)/f0
Ba=-1*gamma*(sin(2.0*X)-2.0*X)/(2.0*X**2.0)
Ecp=Ec*(1-2*Ba/wq0)
Ecvec=-Ecp*(6.0*nvec**2+6.0*nvec+3.0)/12.0
wTvec = -Ej + sqrt(8.0*Ej*Ecp)*(nvec+0.5)+Ecvec
#print diff(wTvec)
if phi==0.4:
print Ba/wq0
if 0:
zr=zeros(N)
X=Np*pi*(diff(wTvec)-f0)/f0
Ba=-0*gamma*(sin(2.0*X)-2.0*X)/(2.0*X**2.0)
ls=-Ba#/(2*Ct)/(2*pi)#/1e6
zr[1:]=ls
#print Ga0/(2*Ct)/(2*pi)
#print zr
wTvec=wTvec+zr
wTvec = wTvec-wdvec
#print wTvec+diff(wTvec)
transmon_levels = Qobj(diag(wTvec[range(N)]))
H=transmon_levels +Omega_vec #- 0.5j*(Omega_true*adag - conj(Omega_true)*a)
final_state = steadystate(H, c_op_list)
return expect( tm_l, final_state)
开发者ID:thomasaref,项目名称:TA_software,代码行数:37,代码来源:qutip_gate_Lamb.py
示例15: H
wVals = np.linspace(-2 * np.pi, 6 * np.pi, 100)
# Solving the dynamics of a spin 1/2 in B0 and Brf in the lab frame
def H(w):
""" Hamiltonian in the lab frame. """
return w0 / 2 * qt.sigmaz() + w1 * np.cos(w) * qt.sigmax()
c1 = np.sqrt(1.0/T1) * qt.destroy(2)
c2 = np.sqrt(1.0/T2) / 2 * qt.destroy(2)
c_op_list = [c1,c2]
wHVals = np.linspace(-2 * np.pi, 6 * np.pi, 51)
MxH = [qt.expect(qt.sigmax(), qt.steadystate(H(w), c_op_list, maxiter=100)) for w in wHVals]
MyH = [qt.expect(qt.sigmay(), qt.steadystate(H(w), c_op_list, maxiter=100)) for w in wHVals]
MzH = [qt.expect(qt.sigmaz(), qt.steadystate(H(w), c_op_list, maxiter=100)) for w in wHVals]
# Plot results of the algebraic formula and the solution given by qutip
plt.ylim([-1.0, 1.0])
plt.plot(wVals, Mx(wVals), label='Mx')
plt.plot(wVals, My(wVals), label='My')
plt.plot(wVals, Mz(wVals), label='Mz')
plt.plot(wHVals, MxH, 'bo', label='MxH')
plt.plot(wHVals, MyH, 'go', label='MyH')
plt.plot(wHVals, MzH, 'ro', label='MzH')
开发者ID:p4t48,项目名称:lab_frame,代码行数:30,代码来源:lab_frame.py
示例16: find_expect
def find_expect(phi=0.1, Omega_vec=3.0, wd=wd, use_pinv=True):
Ej = (phi**2)/(8*Ec) #Ejmax*absolute(cos(pi*phi)) #; % Josephson energy as function of Phi.
wTvec = -Ej + sqrt(8.0*Ej*Ec)*(nvec+0.5)+Ecvec
#wdvec=nvec*sqrt(8.0*Ej*Ec)
#wdvec=nvec*wd
wT = wTvec-wdvec
transmon_levels = Qobj(diag(wT[range(N)]))
H=transmon_levels +Omega_vec #- 0.5j*(Omega_true*adag - conj(Omega_true)*a)
L=liouvillian(H, c_op_list) #ends at same thing as L = H_comm + Lindblad_tm ;
rho_ss = steadystate(L) #same as rho_ss_c but that's in column vector form
tr_mat = tensor([qeye(n) for n in L.dims[0][0]])
#N = prod(L.dims[0][0])
A = L.full()
#D, U= L.eigenstates()
#print D.shape, D
#print diag(D)
b = spre(p).full()-spost(p).full()
a2 = spre(a).full()
tr_vec = transpose(mat2vec(tr_mat.full()))
rho = transpose(mat2vec(rho_ss.full()))
I = identity(N * N)
P = kron(transpose(rho), tr_vec)
Q = I - P
#wp=4.9
#w=-(wp-wd)
spectrum = zeros(len(wlist2), dtype=complex)
for idx, w in enumerate(wlist2):
if use_pinv:
MMR = pinv(-1.0j * w * I + A) #eig(MMR)[0] is equiv to Dint
else:
MMR = dot(Q, solve(-1.0j * w * I + A, Q))
#print diag(1.0/(1j*(wp-wd)*ones(N**2)+D)) #Dint = diag(1./(1i*(wp-wd)*diag(eye(dim^2)) + diag(D)))
#print 1.0/(1j*(wp-wd)*ones(N**2)+D) #Dint = diag(1./(1i*(wp-wd)*diag(eye(dim^2)) + diag(D)))
#U2=squeeze([u.full() for u in U]).transpose()
#Dint=eig(MMR)[0]
#print "MMR", eig(MMR)[1]
#print "Umult", U2*Dint*inv(U2)
s = dot(tr_vec,
dot(a2, dot(MMR, dot(b, transpose(rho)))))
#spectrum[idx] = -2 * real(s[0, 0])
spectrum[idx]=1j*s[0][0] #matches Chi_temp result #(1/theta_steps)
return spectrum
#final_state = steadystate(H, c_op_list)
#print H.shape
#print dir(H)
#U,D = eig(H.full())
#print D
#Uinv = Qobj(inv(D))
#U=Qobj(D)
# Doing the chi integral (gives the susceptibility)
#Dint = 1.0/(1.0j*(wp-wd)) #Qobj(1.0/(1.0j*(wp-wd)*diag(qeye(N**2))))# + diag(D))))
#Hint = H.expm() #U*H*Uinv
#Chi_temp(i,j) += (1.0/theta_steps)*1j*trace(reshape(tm_l*Lint*(p_l - p_r)*rho_ss_c,dim,dim))
#exp1=correlation_2op_1t(H, None, wlist2, c_op_list, a, p, solver="es")
#exp2=correlation_2op_1t(H, None, wlist2, c_op_list, p, a, solver="es", reverse=True)
#exp1=correlation_2op_1t(H, None, wlist2, c_op_list, a, p_l-p_r, solver="es")
#exp1=correlation_ss(H, wlist2, c_op_list, a, p)
#exp2=correlation_ss(H, wlist2, c_op_list, p, a, reverse=True)
exp1=spectrum(H, wlist2, c_op_list, a, p, solver="pi", use_pinv=False)
exp2=spectrum(H, wlist2, c_op_list, p, a, solver="pi", use_pinv=False)
return exp1-exp2
return expect( a, final_state) #tm_l
开发者ID:thomasaref,项目名称:TA_software,代码行数:79,代码来源:qutip_twotone_fixed.py
示例17: conj
print Omega_vec
H=transmon_levels +Omega_vec #- 0.5j*(Omega_true*adag - conj(Omega_true)*a)
print H
#Ht = H.data.T
#from qutip.cy.spmath import zcsr_kron
#from qutip.fastsparse import fast_csr_matrix, fast_identity
#spI=fast_identity(N)
#data = -1j * zcsr_kron(spI, H.data)
#data += 1j * zcsr_kron(Ht, spI)
#print data.toarray()
L=liouvillian(H, c_op_list) #ends at same thing as L = H_comm + Lindblad_tm ;
#print L
#print dir(L)
rho_ss = steadystate(L) #same as rho_ss_c but that's in column vector form
print rho_ss
wlist = linspace(-500.0, 500.0, 201)
use_pinv=True
tr_mat = tensor([qeye(n) for n in L.dims[0][0]])
N = prod(L.dims[0][0])
A = L.full()
D, U= L.eigenstates()
print D.shape, D
print diag(D)
b = spre(p).full()-spost(p).full()
a = spre(a).full()
开发者ID:thomasaref,项目名称:TA_software,代码行数:31,代码来源:qutip_compare.py
示例18: H
# Calculating steady-state solutions using the qutip library
def H(w):
""" Hamiltonian of a spin 1/2 in B0 + Brf using the RWA. """
return (w0 - w)/2 * qt.sigmaz() + w1/2 * qt.sigmax()
# Collsapse operators which should describe the relaxation of the system
c1 = np.sqrt(1.0/T1) * qt.destroy(2)
c2 = np.sqrt(1.0/T2)/2 * qt.destroy(2)
c3 = np.sqrt(1.0/(T1*T2)) * qt.destroy(2)
c4 = np.sqrt(1.0/T2) /2 * qt.create(2)
c5 = np.sqrt(1.0/(T1*T2)) * qt.create(2)
wHVals = np.linspace(-2 * np.pi, 6 * np.pi, 31)
MxH = [qt.expect(qt.sigmax(), qt.steadystate(H(w), [c1,c2,c3], maxiter=10)) for w in wHVals]
MyH = [qt.expect(qt.sigmay(), qt.steadystate(H(w), [c2,c5], maxiter=10)) for w in wHVals]
MzH = [qt.expect(qt.sigmaz(), qt.steadystate(H(w), [c1], maxiter=10)) for w in wHVals]
wVals = np.linspace(-2 * np.pi, 6 * np.pi, 100)
# Plot results of the algebraic formula and the solution given by qutip
plt.ylim([-1.0, 1.0])
plt.plot(wVals, Mx(wVals), label='Mx')
plt.plot(wVals, My(wVals), label='My')
plt.plot(wVals, Mz(wVals), label='Mz')
plt.plot(wHVals, MxH, 'bo', label='MxH')
plt.plot(wHVals, MyH, 'go', label='MyH')
开发者ID:p4t48,项目名称:rotating_frame,代码行数:31,代码来源:rotating_frame.py
示例19: two_tone
def two_tone(vg):
phi, wd=vg
Ej = Ejmax*absolute(cos(pi*phi))#(phi**2)/(8*Ec) # #; % Josephson energy as function of Phi.
wTvec = -Ej + sqrt(8.0*Ej*Ec)*(nvec+0.5)+Ecvec
wdvec=nvec*wd
wT = wTvec-wdvec
transmon_levels = Qobj(diag(wT[range(N)]))
for Om in Omega_op:
H = transmon_levels + Om
exp1=spectrum(H, wlist2, c_op_list, a, p, solver="pi", use_pinv=False)
exp2=spectrum(H, wlist2, c_op_list, p, a, solver="pi", use_pinv=False)
return exp1-exp2
#H_comm = -1j*(kron(It,H) - kron(transpose(H),It));
#L = H_comm + Lindblad_tm + Lindblad_tp #+ Lindblad_deph;
#print L
#L2 = [reshape(eye(N),1,N**2); L]; %Add the condition trace(rho) = 1
#rho_ss = L2\[1;zeros(dim^2,1)]; % Steady state density matrix
#rho_ss_c = rho_ss; %Column vector form
L = liouvillian(H, c_op_list)
#D, V=L.eigenstates()
#print "D", D.shape
#print "V", V.shape
#print L.diag()
#print L.diag().shape
#raise Exception
A = L.full()
rho_ss = steadystate(L)
rho = transpose(mat2vec(rho_ss.full()))
P = kron(transpose(rho), tr_vec)
Q = I - P
if 1:
w=wp-wd
#MMR = pinv(-1.0j * w * I + A)
MMR = dot(Q, solve(-1.0j * w * I + A, Q))
#Lexp=L.expm()
#MMR=Lexp*MMR*Lexp.dag()
# MMR = np.dot(Q, np.linalg.solve(-1.0j * w * I + A, Q))
#print p_l.shape
#print p_l
#print transpose(rho).shape
#print dot(p_l.full(), transpose(rho))
s = dot(tr_vec,
dot(tm_l, dot(MMR, dot(p_l_m_p_r, transpose(rho)))))
return -2 * real(s[0, 0])
D,U = eig(L.full())
Uinv = inv(U)
Dint = diag(1.0/(1.0j*(wp-wd)*Idiag + diag(D)))
Lint = U*Dint*Uinv
#H_comm = -1.0j*(kron(It,H) - kron(transpose(H),It))
#print H
#print H.shape, H.full().shape
#print H_comm, H_comm.shape
#L = H_comm + Lindblad_tm + Lindblad_tp + Lindblad_deph
#p_l=spre(p)
#p_r=spost(p)
#tm_l=spre(a)
#print L.shape, A.shape,
#print dir(L)
#print L.eigenenergies().shape
#print L.eigenstates()#.shape
#L2 = [reshape(eye(dim),1,dim^2);L]; %Add the condition trace(rho) = 1
#rho_ss = L2\[1;zeros(dim^2,1)]; % Steady state density matrix
#rho_ss_c = rho_ss; %Column vector form
#print "tm_l", tm_l.full().shape
#print "Lint", Lint.shape
#print "pl", p_l.full().shape
#print "pr", p_r.full().shape
#print rho.shape, to_super(rho_ss).full().shape
#print "rho", to_super(rho_ss) #rho_ss.full().shape #rho.shape
#print (tm_l.full()*Lint*(p_l.full() - p_r.full())*rho).shape
#raise Exception
return 1j*trace(tm_l*Lint*(p_l_m_p_r)*rho_ss)
开发者ID:thomasaref,项目名称:TA_software,代码行数:97,代码来源:qutip_twotone3.py
注:本文中的qutip.steadystate函数示例由纯净天空整理自Github/MSDocs等源码及文档管理平台,相关代码片段筛选自各路编程大神贡献的开源项目,源码版权归原作者所有,传播和使用请参考对应项目的License;未经允许,请勿转载。 |
请发表评论