I have been coding to calculate the following equation:
Where rho(t) is calculated by solving a quantum mechanics equation. <3| is a row vector with 1 at the 3rd element. |3> is a column vector with 1 at the 3rd element.
To solve the quantum mechanics equation (for the physics experts, I am solving the Lindblad Equation), I am using the qutip module.
The problem is, whenever I do the integration (using the scipy.integrate.quad function), I get a number like -0.842371561579 which is not what I expected (I am hoping for a number like 0.99). Also, I get a warning:
/usr/local/lib/python2.7/dist-packages/scipy/integrate/quadpack.py:352: IntegrationWarning: The integral is probably divergent, or slowly convergent.
warnings.warn(msg, IntegrationWarning)
I tried playing around with the values of the constants in the gamma constant. However, it didn't change any of the output. I have noticed it has to do with calculating rho(t), but that's all I could discover.
What is wrong with the code?
from qutip import *
import numpy as np
import scipy
from scipy.constants import *
hamiltonian = np.array([[215, -104.1, 5.1, -4.3 ,4.7,-15.1 ,-7.8 ],
[-104.1, 220.0, 32.6 ,7.1, 5.4, 8.3, 0.8],
[ 5.1, 32.6, 0.0, -46.8, 1.0 , -8.1, 5.1],
[-4.3, 7.1, -46.8, 125.0, -70.7, -14.7, -61.5],
[ 4.7, 5.4, 1.0, -70.7, 450.0, 89.7, -2.5],
[-15.1, 8.3, -8.1, -14.7, 89.7, 330.0, 32.7],
[-7.8, 0.8, 5.1, -61.5, -2.5, 32.7, 280.0]])
H=Qobj(hamiltonian) #This just converts hamiltonian into the qutip format
ground=H.groundstate()[1]
rho0 = ground*ground.dag()
gamma=(2*pi)*(296*0.695)*(35/150)
print gamma
#collapse operators:
L1 = Qobj(np.array([[1,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L2 = Qobj(np.array([[0,0,0,0,0,0,0],[0,1,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L3 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,1,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L4 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,1,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L5 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,1,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L6 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,1,0],[0,0,0,0,0,0,0]]))
L7 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,1]]))
#Since our gamma variable cannot be directly applied onto the Lindblad operator, we must multiply it with the collapse operators:
L1=gamma*L1
L2=gamma*L2
L3=gamma*L3
L4=gamma*L4
L5=gamma*L5
L6=gamma*L6
L7=gamma*L7
options = Options(nsteps=10000000, atol=1e-20)
#This code is for function-based integration
def integratefunc(x):
results = mesolve(H, rho0, [x], [L1,L2,L3,L4,L5,L6,L7],[], options=options) #this function find rho(t) and evaluates it at t=x
bra3= [[0,0,1,0,0,0,0]] #<3|
bra3q=Qobj(bra3) #convert to qutip format
ket3= [[0],[0],[1],[0],[0],[0],[0]] #|3>
ket3q=Qobj(ket3) #convert to qutip format
tempq=bra3q*(results.states[0]*ket3q) #multiply <3|rho(x)|3>
y=tempq.data.data[0].real #gets the real part of the number, since the imaginary part is zero
return y #returns the value
efficiency, error = scipy.integrate.quad(integratefunc, 0, np.inf) #does the integration to calculate the efficiency and the error
print efficiency, error #prints the final results
