output.states not showing anything for multi-qubit system

[anir...@gmail.com]

I have a code for determining the final state for a 4 qubit system. I am expecting a final state. But, output.states does not show anything. Yet, I am getting graphs properly for occupation prob. of the qubits. Here is my code and error.

 # -*- coding: utf-8 -*-
"""
Created on Thu Aug 24 19:41:32 2017
 
@author: anirban
"""
import matplotlib.pyplot as plt
import numpy as np
import math
from scipy.fftpack import rfft, irfft, fftfreq, fft
import scipy.fftpack
 
from qutip import *
 
def Time_of_swap(wq,w1,w2,w3,g1,g2,g3,wmode,N,psi0,A,tlist,num,T,psi_final):
    
 omega=wmode-wq #the modulation frequency.
 
 # operators
 b1=tensor(qeye(2),destroy(N),qeye(2),qeye(2))
 b2=tensor(qeye(2),qeye(2),destroy(N),qeye(2))
 b3=tensor(qeye(2),qeye(2),qeye(2),destroy(N))
 sm =tensor(destroy(2),qeye(N),qeye(N),qeye(N))
 
 
 # Hamiltonian
 H0=w1*b1.dag()*b1+w2*b2.dag()*b2+w3*b3.dag()*b3+g1*(b1+b1.dag())*(sm+sm.dag())+g2*(b2+b2.dag())*(sm+sm.dag())+g3*(b3+b3.dag())*(sm+sm.dag())
 H1=0.5*tensor(sigmaz(),qeye(2),qeye(2),qeye(2))
 
 c_ops = []
 args1={'omega':omega,'A':A}
 
 #Forming the drive in the transformed basis
 def H1_coeff(t, args1):
  return -wq-(args1['A']*np.sin(args1['omega']*t))
 
 H=[H0,[H1,H1_coeff]]
 
 output = mesolve(H, psi0, tlist, c_ops, [sm.dag() * sm,b1.dag() * b1,b2.dag() *b2,b3.dag()*b3], args1)
 
 print len(output.states)
 
 n_q=output.expect[0]
 n_r1 = output.expect[1]
 n_r2 = output.expect[2]
 n_r3 = output.expect[3]
 
 
 if(wmode==w1):
     n_mode=n_r1
 elif(wmode==w2):
     n_mode=n_r2
 elif(wmode==w3):
     n_mode=n_r3
     
 psi_final=output.states[len(output.states)-1]

         
 fig, axes = plt.subplots(1,1, figsize=(10,6))
 axes.plot(tlist,n_q, label="transmon qubit")
 axes.plot(tlist,n_mode,label="Accessed resonator mode")
 axes.grid()
 #axes.set_xlim(0,30)
 axes.legend(loc=0)
 axes.set_xlabel('Time')
 axes.set_ylabel('Occupation probability')
 axes.set_title('Stimulated Vacuum rabi oscillations')
 plt.show()
 
 
 
#Data
wq = 4.28*2*np.pi  #qubit frequency between |g> and |e>
wq2= 5.0*2*np.pi
w1 = 6.0*2*np.pi   #mode1 frequency
w2 = 6.5*2*np.pi   #mode2 frequency
w3 = 7.0*2*np.pi   #mode3 frequency
g1 = 0.10 * 2 *np.pi  # coupling strength of mode1
g2 = 0.15 * 2 *np.pi  # coupling strength of mode2
g3 = 0.20 * 2 *np.pi  # coupling strength of mode3
 
num=500
T=0.2
#tlist = np.linspace(0,num*T,num) #Duration of pulse
N=2 #no.of fock states
# intial state
psi0 =tensor(basis(2,0),basis(N,1),(basis(N,0)+basis(N,1)).unit(),basis(N,0))
A=5.0*2*np.pi # change this variable

tlist=np.linspace(0,7.11428571429,100)
Time_of_swap(wq,w1,w2,w3,g1,g2,g3,w1,N,psi0,A,tlist,num,T,psi0)

Error:
Traceback (most recent call last):
  File "two_mode.py", line 118, in <module>
    Time_of_swap(wq,w1,w2,w3,g1,g2,g3,w1,N,psi0,A,tlist,num,T,psi0)
  File "two_mode.py", line 41, in Time_of_swap
    print output.states[len(output.states)-1]
IndexError: list index out of range

[Paul Nation]

You are are asking for expectation values in your call to mesolve.  As such, there are no state vectors returned.

-P