Error with state.unit()

[Maxim Moloshenko]

Context: I try to reproduce the results of "Nonequilibrium delocalization-localization transition of photons in circuit quantum elecrodynamics" by Schmidt, Gerace and so on.
Error:
Traceback (most recent call last):
 File "/home/user-admin/Documents/Photons/main.py", line 70, in <module>
state = state.unit()
 File "/usr/local/lib/python2.7/dist-packages/qutip/qobj.py", line 824, in unit
oper_norm=oper_norm, sparse=sparse, tol=tol, maxiter=maxiter)
 File "/usr/local/lib/python2.7/dist-packages/qutip/qobj.py", line 675, in norm
return _sp_L2_norm(self)
 File "/usr/local/lib/python2.7/dist-packages/qutip/sparse.py", line 52, in _sp_L2_norm
return la.norm(op.data.data, 2)
 File "/usr/local/lib/python2.7/dist-packages/scipy/linalg/misc.py", line 121, in norm
return nrm2(a)
_fblas.error: (offx>=0 && offx<len(x)) failed for 2nd keyword offx: dznrm2:offx=0
Script terminated.

I just don't understand how I should tackle the problem. It seems like something deep inside some BLAS  routine, so I am kinda clueless. I would really appreciate if you either told me whats going wrong, or what to do next to find the bug.
Thanks a lot!
Maxim

Here is the essential code:
# Import the amazing QuTiP library
from qutip import *

# Define the start state
state  = (tensor(basis(2,0),fock(20,19),basis(2,0),fock(20,0))).unit()

# Define the Hamilton operator (constant part)
# Local JC Hamiltonian:
omega = 0.4
gg    = 0.6
HJC = omega*tensor(qeye(2),destroy(20)*create(20),qeye(2),qeye(20))\
    + omega*tensor(sigmap()*sigmam(),qeye(20),qeye(2),qeye(20))\
    + gg*(
          tensor(sigmap(),create(20),qeye(2),qeye(20)) +
          tensor(sigmam(),destroy(20),qeye(2),qeye(20))
         )\
    + omega*tensor(qeye(2),qeye(20),qeye(2),destroy(20)*create(20))\
    + omega*tensor(qeye(2),qeye(20),sigmap()*sigmam(),qeye(20))\
    + gg*(
          tensor(qeye(2),qeye(20),sigmap(),create(20)) +
          tensor(qeye(2),qeye(20),sigmam(),destroy(20))
         )

# Coupling part
J = 0.7
temp = tensor(qeye(2),destroy(20),qeye(2),create(20))
Ham = HJC - J*(temp + temp.conj())

# Simulation of the time development, given a certain start state
# For every timestep, calculate exp(-i H t)|state>, and measure the
# number of photons in the left cavity
Hamtt = qeye(1600)
for t in xrange(0,100):
  Hamtt = - 1j * (t) * Ham
  state = Hamtt.expm()*state
  state = state.unit()
  print 'Number of call:',t
  print 'Scalar product:',(state.dag()*state).diag()
  print 'Photons left:',\
        (state.dag()*tensor(qeye(2),num(20),qeye(2),qeye(20))*state).diag()
  print 'Qubit state left:',\
        (state.dag()*tensor(num(2),qeye(20),qeye(2),qeye(20))*state).diag()
  print 'Photons right:',\
        (state.dag()*tensor(qeye(2),qeye(20),qeye(2),num(20))*state).diag()
  print 'Qubit state right:',\
        (state.dag()*tensor(qeye(2),qeye(20),num(2),qeye(20))*state).diag()

[Paul Nation]

Maxim,

Your code seems to work fine on my computer.  Here is a sample of the output:

Number of call: 0
Scalar product: [ 1.]
Photons left: [ 19.]
Qubit state left: [ 0.]
Photons right: [ 0.]
Qubit state right: [ 0.]
Number of call: 1
Scalar product: [ 1.]
Photons left: [ 6.24051947]
Qubit state left: [ 0.66321109]
Photons right: [ 11.30688868]
Qubit state right: [ 0.78938076]
Number of call: 2
Scalar product: [ 1.]
Photons left: [ 1.39624029]
Qubit state left: [ 0.0975253]
Photons right: [ 16.64655971]
Qubit state right: [ 0.85967471]
Number of call: 3
Scalar product: [ 1.]
Photons left: [ 0.37475689]
Qubit state left: [ 0.25045689]
Photons right: [ 17.77055895]
Qubit state right: [ 0.60422727]

It does appear that you have a problem in your BLAS implementation.  In order to fix it, I would think you need to reinstall your BLAS libraries.  Although it could also be an issue with NumPy passing an illegal value to the dznrm2 function from numpy.norm.

Also, just a tip, instead of always writing create(20), destroy(20),….  why not just define them once like

a=destroy(20)

ad=create(20)  [or ad.a.dag()]

and then you will save your self lots of time, and your code will look more physical.  Also, it is best to stick to using only basis or fock and not mix the two.  They do the exact same thing, but it makes reading your code easier.

Paul

---------------------------------------
Paul D. Nation
Assistant Professor
Korea University
Department of Physics
Anam-dong 5, Seongbuk-gu
Seoul 136-713, South Korea
Email: pna...@korea.ac.kr
Phone: +82-02-3290-3092
Web:   nqdl.korea.ac.kr
---------------------------------------

--
You received this message because you are subscribed to the Google Groups "qutip" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qutip+un...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[jrjoh...@gmail.com]

Hi Maxim

I could reproduce your problem on my computer, and it occurs after 10 or so iteration of the loop. The problem is that the state vector in your calculation end up having no nonzero entries (which in turn probably is due to that you take too large time steps, so that the state vector normalization fails). For example, after 10 iterations I get

In [8]: state.data.getnnz()

Out[8]: 0

In [9]: state.data.data

Out[9]: array([], dtype=complex128)

In [11]: import scipy.linalg as la

            la.norm(state.data.data)

---------------------------------------------------------------------------

error                                     Traceback (most recent call last)

<ipython-input-11-21e2d3203a6d> in <module>()

      1 import scipy.linalg as la

      2 

----> 3 la.norm(state.data.data)

/usr/local/lib/python3.3/dist-packages/scipy/linalg/misc.py in norm(a, ord)

    120         func_name = _nrm2_prefix.get(a.dtype.char, 'd') + 'nrm2'

    121         nrm2 = getattr(blas, func_name)

--> 122         return nrm2(a)

    123     return np.linalg.norm(a, ord=ord)

    124 

error: (offx>=0 && offx<len(x)) failed for 2nd keyword offx: dznrm2:offx=0

I guess the way this error is presented could be considered a bug in scipy (empty array input to scipy.linalg.norm should probably give a nicer error message) but the problem with your programs is that the state vector becomes an empty array in the first place. 

To avoid that I think you'd need to take much smaller time steps in the evolution. Or why not use for example mesolve or mcsolve instead of explicitly exponentiating the Hamiltonian. That way you can also include dissipation in the problem.

In fact, a couple of months back I also had a quick look at reproducing some of the results from the paper you are studying. I didn't pursue it very far, but I've put up the notebook I worked on as a gist:

https://gist.github.com/jrjohansson/5963035

http://nbviewer.ipython.org/5963035

Since you are looking at that paper now perhaps this calculation could be of some interest to you. At least it demonstrates how to set up the problem studied in Schmidt et al, in qutip. It takes a quite long time to run so have patience..

Best regards

Rob