The use of collapse operator for solving time-evolution operator

[j794...@uwaterloo.ca]

I am trying to calculate time-evolution operator of some processes (like the evolution of JC model). If we want to add dissipation into my system, typically, the dynamics can be well solved by 'mesolve'. 

But now I would like to get the time-evolution operator so I use 'propagator' function and input some collapse operators like c_op=[a1, a2...]. It seems that the c_ops cannot affect the result in any way, for the function gives me the same U matrix with any different c_ops. 

Is this a QuTip bug? Or it is just because I did not apply the correct solver Options?

[j794...@uwaterloo.ca]

Dear all,

Have anyone met this problem?

[Andrew M. C. Dawes]

It’s hard to help without at least an example section of the code you are trying to run.

If you can post something with and without c_ops that returns the same result, it is more likely someone can identify the cause.

--
You received this message because you are subscribed to the Google Groups "QuTiP: Quantum Toolbox in Python" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qutip+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/qutip/3f3d8629-1694-49e5-8406-06319186db21n%40googlegroups.com.

[j794...@uwaterloo.ca]

I tried to write an example to explain my question, but another error shows up. Here is the code and error

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

wc = 1.0 * 2 * np.pi

wa = 1.0 * 2 * np.pi

g = 0.05 * 2 * np.pi

kappa = 0.005

gamma = 0.05

N = 15

step = 1001

vac = tensor(basis(N,0),basis(2,0))

a = tensor(destroy(N),qeye(2))

sm = tensor(qeye(N),destroy(2))

H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())

c_ops = [np.sqrt(kappa)*a,np.sqrt(gamma)*sm]

t_list = np.linspace(0,25.0,step)

psi0 = sm.dag() * vac

result1 = mesolve(H,psi0,t_list,c_ops,[a.dag()*a,sm.dag()*sm])

result0 = mcsolve(H,psi0,t_list,[],[a.dag()*a,sm.dag()*sm])

n_c_0 = result0.expect[0]

n_a_0 = result0.expect[1] 

n_c_1 = result1.expect[0]

n_a_1 = result1.expect[1]

fig,axes = plt.subplots(1,1,figsize = (10,6))

axes.plot(t_list,n_c_0,label = 'Cavity(without dissipasion)')

axes.plot(t_list,n_a_0,label = 'Atom(without dissipasion)')

axes.plot(t_list,n_c_1,'-.',label = 'Cavity(with dissipasion)')

axes.plot(t_list,n_a_1,'-.',label = 'Atom(with dissipasion)')

axes.legend(loc = 'best')

axes.set_xlabel('Time')

axes.set_ylabel('Occupation probability')

axes.set_title('Vacuum Rabi oscillations')

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

This part works well. But when I tried to solve the time-evolution operator:

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

U1 = propagator(H, t_list,c_ops)

U0 = propagator(H, t_list, [])

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-4-8134c024ad2f> in <module> ----> 1 U1 = propagator(H, t_list,c_ops) 2 U0 = propagator(H, t_list, []) C:\anaconda3\lib\site-packages\qutip\propagator.py in propagator(H, t, c_op_list, args, options, unitary_mode, parallel, progress_bar, _safe_mode, **kwargs) 239 rho0 = Qobj(sp.csr_matrix(([1], ([row_idx], [col_idx])), 240 shape=(N, N), dtype=complex)) --> 241 output = mesolve(H, rho0, tlist, c_op_list, [], args, options, 242 _safe_mode=False) 243 for k, t in enumerate(tlist): C:\anaconda3\lib\site-packages\qutip\mesolve.py in mesolve(H, rho0, tlist, c_ops, e_ops, args, options, progress_bar, _safe_mode) 263 raise Exception("Invalid H type") 264 --> 265 func, ode_args = ss.makefunc(ss, rho0, args, e_ops, options) 266 267 if _safe_mode: C:\anaconda3\lib\site-packages\qutip\mesolve.py in _qobjevo_set(HS, rho0, args, e_ops, opt) 349 # Should be caught earlier in mesolve. 350 raise ValueError("rho0 must be a ket, density matrix or superoperator") --> 351 _test_liouvillian_dimensions(H_td.cte.dims, rho0.dims) 352 return func, () 353 C:\anaconda3\lib\site-packages\qutip\mesolve.py in _test_liouvillian_dimensions(L_dims, rho_dims) 330 raise ValueError("Liouvillian had nonsquare dims: " + str(L_dims)) 331 if not ((L_dims[1] == rho_dims) or (L_dims[1] == rho_dims[0])): --> 332 raise ValueError("".join([ 333 "incompatible Liouvillian and state dimensions: ", 334 str(L_dims), " and ", str(rho_dims), ValueError: incompatible Liouvillian and state dimensions: [[[15, 2], [15, 2]], [[15, 2], [15, 2]]] and [[30], [30]]

[Andrew M. C. Dawes]

Hmm, I've run the first cell (success), and then immediately after, the two U1 and U0 lines without any error. I'm using the following:

Software	Version
QuTiP	4.5.3
Numpy	1.20.2
SciPy	1.6.2
matplotlib	3.3.4
Cython	0.29.23
Number of CPUs	2
BLAS Info	INTEL MKL
IPython	7.22.0
Python	3.7.10 (default, Feb 26 2021, 10:16:00) [Clang 10.0.0 ]
OS	posix [darwin]
Wed Nov 10 06:52:21 2021 PST

To view this discussion on the web visit https://groups.google.com/d/msgid/qutip/87a04031-6a06-470d-833b-61e389c43826n%40googlegroups.com.

[j794...@uwaterloo.ca]

Thank you for checking the code!

Also, my setup is:

Software               Version
QuTiP                        4.6.2
Numpy                     1.20.1

SciPy                        1.6.2
matplotlib                3.3.4
Cython                    0.29.23

Number of CPUs         6

BLAS Info          INTEL MKL
IPython                   7.22.0

Python               3.8.8 (default, Apr 13 2021, 15:08:03) [MSC v.1916 64 bit (AMD64)]
OS                          nt [win32]

It is really weird that latest version performs worse than previous version. By the way, what is the output of U1 and U0 in your results?

U0[-1] and U1[-1] make more sense though. 

[j794...@uwaterloo.ca]

Maybe I found where the incompatibility comes from. In this part of the propagator function:

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

        else:

            progress_bar.start(N * N)

            for n in range(N * N):

                progress_bar.update(n)

                col_idx, row_idx = np.unravel_index(n, (N, N))

                rho0 = Qobj(sp.csr_matrix(([1], ([row_idx], [col_idx])),

                                        shape=(N, N), dtype=complex))

                output = mesolve(H, rho0, tlist, c_op_list, [], args, options,

                                 _safe_mode=False)

                for k, t in enumerate(tlist):

                    u[:, n, k] = mat2vec(output.states[k].full()).T

            progress_bar.finished()

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

The dims of rho0 was not correctly initialized. It was not always compatible with H when we study multi-body system. And when I tried to solve it like:

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

else:

            progress_bar.start(N * N)

            for n in range(N * N):

                progress_bar.update(n)

                col_idx, row_idx = np.unravel_index(n, (N, N))

                rho0 = Qobj(sp.csr_matrix(([1], ([row_idx], [col_idx])),

                                        shape=(N, N), dtype=complex), dims=dims)

                output = mesolve(H, rho0, tlist, c_op_list, [], args, options,

                                 _safe_mode=False)

                for k, t in enumerate(tlist):

                    u[:, n, k] = mat2vec(output.states[k].full()).T

            progress_bar.finished()

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

or like:

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

else:

            progress_bar.start(N * N)

            for n in range(N * N):

                progress_bar.update(n)

                col_idx, row_idx = np.unravel_index(n, (N, N))

                rho0 = Qobj(sp.csr_matrix(([1], ([row_idx], [col_idx])),

                                        shape=(N, N), dtype=complex))

rho0.dims = dims

                output = mesolve(H, rho0, tlist, c_op_list, [], args, options,

                                 _safe_mode=False)

                for k, t in enumerate(tlist):

                    u[:, n, k] = mat2vec(output.states[k].full()).T

            progress_bar.finished()

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

After editing the .py file in library and run my code, it always tells me the notebook ends.

[j794...@uwaterloo.ca]

The accuracy of propagator may not be good. I compared the result following the above codes like:

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

U0 = propagator(H, t_list, [], unitary_mode='batch')

U0[-1]*psi0

U0_ = propagator(H, t_list, [], unitary_mode='single')

U0_[-1]*psi0

result0.states[-1]

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

These three results are totally different.  

[Andrew M. C. Dawes]

You probably don't want to be forcing dimensions inside that function, that would create some issues that would be really hard to track down later on. I suspect something is off about how the dimensions are being defined before it goes to that function. I'll take another look after upgrading to 4.6.2 to see if I can replicate it after that.

To view this discussion on the web visit https://groups.google.com/d/msgid/qutip/e294d403-df57-4f61-a841-655efa4610e5n%40googlegroups.com.

[Andrew M. C. Dawes]

Ok, I see what you mean. This is a new rho0 that is created inside the function, but was being created with incompatible dimensions. That sounds like a good candidate for an issue in github. I'm happy to submit a pull request to see about fixing this, it does look like an easy fix, but a PR is the way to get it seen and fixed officially, let me know if you would prefer to do so.

[Andrew M. C. Dawes]

Ok, the dims issue is already fixed in the code on github, so no need to file a bug. Next release should include it.

For reference, here is the issue and discussion. You may find the comments about numerical tolerances useful... see if that helps with your observation of the difference in the results you mentioned above.

https://github.com/qutip/qutip/issues/1585

[j794...@uwaterloo.ca]

That really helps. I tried the four ways of calculating final state(mcsolve, mesolve, propagator single, propagator batch). Even though the results are still not completely equal, but taking the number to 3 decimals is now doable. 

[j794...@uwaterloo.ca]

The calculation behind the mesolve and propagator is the key conception of qutip, as my personal perspective. It is basically the conversations between Hilbert and Liouville. I am still trying to understand this interesting structure transform in more depth. 

Back to the dimension compatibility, I found developers have changed the code in propagator.py for 4.7 version. I will just download it and have the one in 4.6 replaced. 

Then maybe we can talk about my original question, the collapse is not working in propagator. 

[Simon Cross]

Hi Chen and Andrew,

Andrew, thank you for digging into this and finding the existing fix -- https://github.com/qutip/qutip/pull/1588.

Chen -- sorry you tripped over the bug again. I'm hoping that 4.6.3 will be released later this month.

Chen: Regarding your original example. You said that U0 and U1 were the same, but when I run your example against the fixed QuTiP, the elements of U0 and U1 aren't even the same kind of Qobj. The U0 elements are operators and the U1 elements are superoperators as expected. If I adjust the collapse operators, the items of U1 change as expected. Are you still experiencing an issue?

Yours sincerely,

Simon Cross

[j794...@uwaterloo.ca]

I have also tried my example, and that works!

But sadly in my project code, it still does not work. The bug is just as I described at the top. I will try to figure this out. If there is something interesting point,  I will share it here.

Thank you all for your support! 

[Simon Cross]

Pleasure!

[j794...@uwaterloo.ca]

Sadly, the propagator is still not working in my project code.

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

H = H_data_generator(t_list2, pulse_paras2, device_paras2)

E1 = propagator(H, t_list2, progress_bar=True, c_ops_list=[], \

                                        options=Options(rhs_with_state=True))

E2 = propagator(H, t_list2, progress_bar=True, c_ops_list=c_ops_list, \

                                        options=Options(rhs_with_state=True))

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

And E1==E2 returns True for all the elements of the array. 

The H is generated by my self-defined function, which act as a pulse generator to change my time-dependent Hamiltonian. The 'cubic' interpolation was applied in the data, by the way. I think these are just typical QuTip style codes. I am really confused about it. 

I will try different Hamiltonian and see what happens. 

在2021年11月16日星期二 UTC+8 下午5:28:06<hodg...@gmail.com> 写道：

Pleasure!

[j794...@uwaterloo.ca]

The problem is on the flow control.

the program goes to this stream：

if len(c_op_list) == 0 and H0.isoper:

even a none [ ] collapse list was given.

I am trying to figure out why this happens.

[j794...@uwaterloo.ca]

I found the issue. 

The parameter for collapse operators in propagator is 'c_op_list'  but not 'c_ops' like it in mesolve.

I always input it by 'c_ops_list', following the style of mesolve. 

Because the '**arwgs' exists, a unexpected parameter input will not raise error.

[Andrew M. C. Dawes]

Ahh, good catch, yes that difference in variable names is not ideal. In a quick search of the repo, there are some functions that use c_ops=[] and some that use c_op_list=[] (and even one that does use c_ops_list=[] (piqs.py).

I'm glad you got it sorted out, this may be a good low-priority issue as it would be nice to standardize these naming conventions.

--

You received this message because you are subscribed to the Google Groups "QuTiP: Quantum Toolbox in Python" group.
To unsubscribe from this group and stop receiving emails from it, send an email to qutip+un...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/qutip/d147f64a-7493-4908-bc1b-3daf59661c63n%40googlegroups.com.

[Asier Galicia]

Hey, I posted an issue in the QuTiP repository summarizing what was discussed here. I hope it is useful!

Best regards,

Asier Galicia.