Python qutip.mesolve函数代码示例

def test_ssesolve_homodyne():
    "Stochastic: smesolve: homodyne"
    tol = 0.01

    N = 4
    gamma = 0.25
    ntraj = 25
    nsubsteps = 100
    a = destroy(N)

    H = a.dag() * a
    psi0 = coherent(N, 0.5)
    sc_ops = [np.sqrt(gamma) * a]
    e_ops = [a.dag() * a, a + a.dag(), (-1j) * (a - a.dag())]

    times = np.linspace(0, 2.5, 50)
    res_ref = mesolve(H, psi0, times, sc_ops, e_ops)
    res = smesolve(
        H,
        psi0,
        times,
        [],
        sc_ops,
        e_ops,
        ntraj=ntraj,
        nsubsteps=nsubsteps,
        method="homodyne",
        store_measurement=True,
        map_func=parallel_map,
    )

    assert_(all([np.mean(abs(res.expect[idx] - res_ref.expect[idx])) < tol for idx in range(len(e_ops))]))

    assert_(len(res.measurement) == ntraj)
    assert_(all([m.shape == (len(times), len(sc_ops)) for m in res.measurement]))

def test_ssesolve_heterodyne():
    "Stochastic: ssesolve: heterodyne, time-dependent H"
    tol = 0.01

    N = 4
    gamma = 0.25
    ntraj = 25
    nsubsteps = 100
    a = destroy(N)

    H = [[a.dag() * a,f]]
    psi0 = coherent(N, 0.5)
    sc_ops = [np.sqrt(gamma) * a, np.sqrt(gamma) * a*0.5]
    e_ops = [a.dag() * a, a + a.dag(), (-1j)*(a - a.dag())]

    times = np.linspace(0, 2.5, 50)
    res_ref = mesolve(H, psi0, times, sc_ops, e_ops, args={"a":2})
    res = ssesolve(H, psi0, times, sc_ops, e_ops,
                   ntraj=ntraj, nsubsteps=nsubsteps,
                   method='heterodyne', store_measurement=True,
                   map_func=parallel_map, args={"a":2})

    assert_(all([np.mean(abs(res.expect[idx] - res_ref.expect[idx])) < tol
                 for idx in range(len(e_ops))]))

    assert_(len(res.measurement) == ntraj)
    assert_(all([m.shape == (len(times), len(sc_ops), 2)
                 for m in res.measurement]))

def carrier_flop(rho0, W, eta, delta, theta, phi, c_op_list = [], return_op_list = []):
    ''' Return values of atom and motion populations during carrier Rabi flop 
    for rotation angles theta. Calls numerical solution of master equation.
    @ var rho0: initial density matrix
    @ var W: bare Rabi frequency
    @ var eta: Lamb-Dicke parameter
    @ var delta: detuning between atom and motion
    @ var theta: list of Rabi rotation angles (i.e. theta, or g*time)
    @ var phi: phase of the input laser pulse
    @ var c_op_list: list of collapse operators for the master equation treatment
    @ var return_op_list: list of population operators the values of which will be returned 
    
    returns: time, populations of motional mode and atom
    '''    
    N = shape(rho0.data)[0]/2 # assume N Fock states and two atom states
    a = tensor(destroy(N), qeye(2))
    Wc = qeye(N)
    Wc.data = csr_matrix( qeye(N).data.dot( np.diag(rabi_coupling(N,0,eta) ) ) )    
    sm = tensor( Wc, destroy(2))

    # use the rotating wave approxiation
    H = delta * a.dag() * a + \
         (1./2.)* W * (sm.dag()*exp(1j*phi) + sm*exp(-1j*phi))
    
    if hasattr(theta, '__len__'):
        if len(theta)>1: # I need to be able to pass a list of length zero and not get an error
            time = theta/W
    else:
            time = theta/W        

    output = mesolve(H, rho0, time, c_op_list, return_op_list)
    return time, output

    def testFloquetUnitary(self):
        """
        Floquet: test unitary evolution of time-dependent two-level system
        """

        delta = 1.0 * 2 * np.pi
        eps0 = 1.0 * 2 * np.pi
        A = 0.5 * 2 * np.pi
        omega = np.sqrt(delta ** 2 + eps0 ** 2)
        T = (2 * np.pi) / omega
        tlist = np.linspace(0.0, 2 * T, 101)
        psi0 = rand_ket(2)
        H0 = - eps0 / 2.0 * sigmaz() - delta / 2.0 * sigmax()
        H1 = A / 2.0 * sigmax()
        args = {'w': omega}
        H = [H0, [H1, lambda t, args: np.sin(args['w'] * t)]]
        e_ops = [num(2)]

        # solve schrodinger equation with floquet solver
        sol = fsesolve(H, psi0, tlist, e_ops, T, args)

        # compare with results from standard schrodinger equation
        sol_ref = mesolve(H, psi0, tlist, [], e_ops, args)

        assert_(max(abs(sol.expect[0] - sol_ref.expect[0])) < 1e-4)

def testHOZeroTemperature():
    """
    brmesolve: harmonic oscillator, zero temperature
    """

    N = 10
    w0 = 1.0 * 2 * np.pi
    g = 0.05 * w0
    kappa = 0.15

    times = np.linspace(0, 25, 1000)
    a = destroy(N)
    H = w0 * a.dag() * a + g * (a + a.dag())
    psi0 = ket2dm((basis(N, 4) + basis(N, 2) + basis(N, 0)).unit())

    c_ops = [np.sqrt(kappa) * a]
    a_ops = [a + a.dag()]
    e_ops = [a.dag() * a, a + a.dag()]

    res_me = mesolve(H, psi0, times, c_ops, e_ops)
    res_brme = brmesolve(H, psi0, times, a_ops, e_ops,
                         spectra_cb=[lambda w: kappa * (w >= 0)])

    for idx, e in enumerate(e_ops):
        diff = abs(res_me.expect[idx] - res_brme.expect[idx]).max()
        assert_(diff < 1e-2)

def rsb_flop(rho0, W, eta, delta, theta, phi, c_op_list = [], return_op_list = []):
    ''' Return values of atom and motion populations during red sideband Rabi flop 
    for rotation angles theta. Calls numerical solution of master equation for the 
    Jaynes-Cummings Hamiltonian.
    @ var rho0: initial density matrix
    @ var W: bare Rabi frequency
    @ var delta: detuning between atom and motion
    @ var theta: list of Rabi rotation angle (i.e. theta, or g*time)
    @ var phi: phase of the input laser pulse
    @ var c_op_list: list of collapse operators for the master equation treatment
    @ var return_op_list: list of population operators the values of which will be returned 
    
    returns: time, populations of motional mode and atom
    '''
    N = shape(rho0.data)[0]/2 # assume N Fock states and two atom states
    a = tensor(destroy(N), qeye(2))
    sm = tensor( qeye(N), destroy(2))
    Wrsb = destroy(N)
    one_then_zero = ([float(x<1) for x in range(N)])
    Wrsb.data = csr_matrix( destroy(N).data.dot( np.diag( rabi_coupling(N,-1,eta) / np.sqrt(one_then_zero+np.linspace(0,N-1,N)) ) ) ) 
    Arsb = tensor(Wrsb, qeye(2))
    # use the rotating wave approxiation
    # Note that the regular a, a.dag() is used for the time evolution of the oscillator
    # Arsb is the destruction operator including the state dependent coupling strength
    H = delta * a.dag() * a + \
        (1./2.) * W * (Arsb.dag() * sm * exp(1j*phi) + Arsb * sm.dag() * exp(-1j*phi))
    
    if hasattr(theta, '__len__'):
        if len(theta)>1: # I need to be able to pass a list of length zero and not get an error
            time = theta/(eta*W)
    else:
            time = theta/(eta*W)        

    output = mesolve(H, rho0, time, c_op_list, return_op_list)
    return time, output

def testHOFiniteTemperatureStates():
    """
    brmesolve: harmonic oscillator, finite temperature, states
    """

    N = 10
    w0 = 1.0 * 2 * np.pi
    g = 0.05 * w0
    kappa = 0.25
    times = np.linspace(0, 25, 1000)
    a = destroy(N)
    H = w0 * a.dag() * a + g * (a + a.dag())
    psi0 = ket2dm((basis(N, 4) + basis(N, 2) + basis(N, 0)).unit())

    n_th = 1.5
    w_th = w0/np.log(1 + 1/n_th)

    def S_w(w):
        if w >= 0:
            return (n_th + 1) * kappa
        else:
            return (n_th + 1) * kappa * np.exp(w / w_th)

    c_ops = [np.sqrt(kappa * (n_th + 1)) * a, np.sqrt(kappa * n_th) * a.dag()]
    a_ops = [a + a.dag()]
    e_ops = []

    res_me = mesolve(H, psi0, times, c_ops, e_ops)
    res_brme = brmesolve(H, psi0, times, a_ops, e_ops, [S_w])

    n_me = expect(a.dag() * a, res_me.states)
    n_brme = expect(a.dag() * a, res_brme.states)

    diff = abs(n_me - n_brme).max()
    assert_(diff < 1e-2)

    def testHOFiniteTemperature(self):
        "brmesolve: harmonic oscillator, finite temperature"

        N = 10
        w0 = 1.0 * 2 * np.pi
        g = 0.05 * w0
        kappa = 0.15
        times = np.linspace(0, 25, 1000)
        a = destroy(N)
        H = w0 * a.dag() * a + g * (a + a.dag())
        psi0 = ket2dm((basis(N, 4) + basis(N, 2) + basis(N, 0)).unit())

        n_th = 1.5
        w_th = w0/np.log(1 + 1/n_th)

        def S_w(w):
            if w >= 0:
                return (n_th + 1) * kappa
            else:
                return (n_th + 1) * kappa * np.exp(w / w_th)

        c_ops = [np.sqrt(kappa * (n_th + 1)) * a, np.sqrt(kappa * n_th) * a.dag()]
        a_ops = [a + a.dag()]
        e_ops = [a.dag() * a, a + a.dag()]

        res_me = mesolve(H, psi0, times, c_ops, e_ops)
        res_brme = brmesolve(H, psi0, times, a_ops, e_ops, [S_w])

        for idx, e in enumerate(e_ops):
            diff = abs(res_me.expect[idx] - res_brme.expect[idx]).max()
            assert_(diff < 1e-2)

def test_fn_list_td_corr():
    """
    correlation: comparing TLS emission correlations (fn-list td format)
    """

    # calculate emission zero-delay second order correlation, g2(0), for TLS
    # with following parameters:
    #   gamma = 1, omega = 2, tp = 0.5
    # Then: g2(0)~0.57
    sm = destroy(2)
    args = {"t_off": 1, "tp": 0.5}
    H = [[2 * (sm+sm.dag()),
          lambda t, args: np.exp(-(t-args["t_off"])**2 / (2*args["tp"]**2))]]
    tlist = linspace(0, 5, 50)
    corr = correlation_3op_2t(H, fock(2, 0), tlist, tlist, [sm],
                              sm.dag(), sm.dag() * sm, sm, args=args)
    # integrate w/ 2D trapezoidal rule
    dt = (tlist[-1]-tlist[0]) / (np.shape(tlist)[0]-1)
    s1 = corr[0, 0] + corr[-1, 0] + corr[0, -1] + corr[-1, -1]
    s2 = sum(corr[1:-1, 0]) + sum(corr[1:-1, -1]) + \
        sum(corr[0, 1:-1]) + sum(corr[-1, 1:-1])
    s3 = sum(corr[1:-1, 1:-1])

    exp_n_in = np.trapz(
        mesolve(
            H, fock(2, 0), tlist, [sm], [sm.dag()*sm], args=args
        ).expect[0], tlist
    )
    # factor of 2 from negative time correlations
    g20 = abs(
        sum(0.5*dt**2*(s1 + 2*s2 + 4*s3)) / exp_n_in**2
    )

    assert_(abs(g20-0.57) < 1e-1)

def test_smesolve_photocurrent():
    "Stochastic: photocurrent_mesolve"
    tol = 0.01

    N = 4
    gamma = 0.25
    ntraj = 20
    nsubsteps = 100
    a = destroy(N)

    H = [[a.dag() * a,f]]
    psi0 = coherent(N, 0.5)
    sc_ops = [np.sqrt(gamma) * a, np.sqrt(gamma) * a * 0.5]
    e_ops = [a.dag() * a, a + a.dag(), (-1j)*(a - a.dag())]

    times = np.linspace(0, 1.0, 21)
    res_ref = mesolve(H, psi0, times, sc_ops, e_ops, args={"a":2})
    res = photocurrent_mesolve(H, psi0, times, [], sc_ops, e_ops, args={"a":2},
                   ntraj=ntraj, nsubsteps=nsubsteps,  store_measurement=True,
                   map_func=parallel_map)

    assert_(all([np.mean(abs(res.expect[idx] - res_ref.expect[idx])) < tol
                 for idx in range(len(e_ops))]))
    assert_(len(res.measurement) == ntraj)
    assert_(all([m.shape == (len(times), len(sc_ops))
                 for m in res.measurement]))

    def simulate(self):
        
        # setup time-dependent Hamiltonians sequentially
        psi0 = self.initial_state
        
        self.sim_states = []
        self.sim_tlist = []
        
        for [H_td, td_coeff, args, tlist] in self._H_td:
            if H_td is not None:
                H = [self.H_ti, [H_td, td_coeff]]
            else:
                H = self.H_ti
                args = None
                
            output = qt.mesolve(H, psi0, tlist, self.c_ops, [], args=args)
            
            # append all but the last step--will the initial step in the next simulation
            self.sim_states.extend(output.states[:-1])
            self.sim_tlist.extend(tlist[:-1])
            
            psi0 = output.states[-1]

        self.sim_states.extend([output.states[-1]])
        self.sim_tlist.extend([tlist[-1]])
            
        self.sim_states = np.array(self.sim_states)
        self.sim_tlist = np.array(self.sim_tlist)

def _qubit_integrate(tlist, psi0, epsilon, delta, g1, g2, solver):

    H = epsilon / 2.0 * sigmaz() + delta / 2.0 * sigmax()

    c_op_list = []

    rate = g1
    if rate > 0.0:
        c_op_list.append(np.sqrt(rate) * sigmam())

    rate = g2
    if rate > 0.0:
        c_op_list.append(np.sqrt(rate) * sigmaz())

    e_ops = [sigmax(), sigmay(), sigmaz()]

    if solver == "me":
        output = mesolve(H, psi0, tlist, c_op_list, e_ops)
    elif solver == "es":
        output = essolve(H, psi0, tlist, c_op_list, e_ops)
    elif solver == "mc":
        output = mcsolve(H, psi0, tlist, c_op_list, e_ops, ntraj=750)
    else:
        raise ValueError("unknown solver")

    return output.expect[0], output.expect[1], output.expect[2]

def testJCZeroTemperature():
    """
    brmesolve: Jaynes-Cummings model, zero temperature
    """

    N = 10
    a = tensor(destroy(N), identity(2))
    sm = tensor(identity(N), destroy(2))
    psi0 = ket2dm(tensor(basis(N, 1), basis(2, 0)))
    a_ops = [(a + a.dag())]
    e_ops = [a.dag() * a, sm.dag() * sm]

    w0 = 1.0 * 2 * np.pi
    g = 0.05 * 2 * np.pi
    kappa = 0.05
    times = np.linspace(0, 2 * 2 * np.pi / g, 1000)

    c_ops = [np.sqrt(kappa) * a]
    H = w0 * a.dag() * a + w0 * sm.dag() * sm + \
        g * (a + a.dag()) * (sm + sm.dag())

    res_me = mesolve(H, psi0, times, c_ops, e_ops)
    res_brme = brmesolve(H, psi0, times, a_ops, e_ops,
                         spectra_cb=[lambda w: kappa * (w >= 0)])

    for idx, e in enumerate(e_ops):
        diff = abs(res_me.expect[idx] - res_brme.expect[idx]).max()
        assert_(diff < 5e-2)  # accept 5% error

    def testMETDDecayAsArray(self):
        "mesolve: time-dependence as array with super as init cond"
        me_error = 1e-5

        N = 10
        a = destroy(N)
        H = a.dag() * a
        psi0 = basis(N, 9)
        rho0vec = operator_to_vector(psi0 * psi0.dag())
        E0 = sprepost(qeye(N), qeye(N))
        kappa = 0.2
        tlist = np.linspace(0, 10, 1000)
        c_op_list = [[a, np.sqrt(kappa * np.exp(-tlist))]]
        out1 = mesolve(H, psi0, tlist, c_op_list, [])
        out2 = mesolve(H, E0, tlist, c_op_list, [])
        fid = self.fidelitycheck(out1, out2, rho0vec)
        assert_(max(abs(1.0 - fid)) < me_error, True)

 def mesolve(self, exps=[]):
     """solve
     Interface to qutip mesolve for the system.
     :param exps: List of expectation values to calculate at
     each timestep.
     Defaults to empty.
     """
     return qt.mesolve(self.hamiltonian()[0], self.initial_state, self.tlist, self._c_ops(), exps)