Python qutip.steadystate函数代码示例

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

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)

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)

    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()

    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

    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

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)

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

    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

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)

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)

    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)

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

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)

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')