Hi,
I've been trying to find the propagator for a time-dependent hamiltonian (with time-dependent collapse operators as well), and have been getting an error saying that the t_list lengths do not match in Qutip 5.0. I also tried to run it without using the collapse operators, but I get the same issue. The same code runs in Qutip 4.7.
My goal is to actually find two-time correlation functions of a time-dependent operator (following this paper:https://journals.aps.org/pra/abstract/10.1103/PhysRevA.102.023717). What I did do (and what has worked) is decompose the operator into time-independent parts and then use Qutip's correlation_2op_2t. This takes a really long time, since I need to do this for 8 different pairs of operators. I thought of instead finding the propagator and manually finding the correlation using it instead.
From the code below, you can see that mesolve works, but the propagator doesn't.
Code to Reproduce the Buggamma_1d_1 = Gamma
gamma_1d_2 = Gamma
## Input pulse
u = input_pulse
mod_u_2 = np.abs(u)**2
integral_u_2 = np.cumsum(mod_u_2) * t[1]
g_u = np.conjugate(u)/np.sqrt(1 - integral_u_2)
## Hamiltonian
emitter_dims = 3
mode_dims = 7
a = 1.1
psi_coherent = coherent(mode_dims, a)
## Mode dims is defined above
w0_1 = 2 * np.pi * 0
alpha_1 = 2 * np.pi * 0.3
w0_2 = 2 * np.pi * 0
alpha_2 = 2 * np.pi * 0.3
gamma_intrinsic = 2 * np.pi * 0.0
sig_1 = tensor(destroy(emitter_dims), qeye(emitter_dims), qeye(mode_dims))
sig_2 = tensor(qeye(emitter_dims), destroy(emitter_dims), qeye(mode_dims))
a_u = tensor(qeye(emitter_dims), qeye(emitter_dims), destroy(mode_dims))
H0 = (w0_1 * sig_1.dag() * sig_1 - alpha_1/2 * sig_1.dag() * sig_1.dag() * sig_1 * sig_1) + \
(w0_2 * sig_2.dag() * sig_2 - alpha_2/2 * sig_2.dag() * sig_2.dag() * sig_2 * sig_2) + \
0.5j * np.sqrt(gamma_1d_1 * gamma_1d_2) * (sig_1.dag() * sig_2 - sig_2.dag() * sig_1)
H_int1 = [0.5j * (a_u.dag() * (np.sqrt(gamma_1d_1) * sig_1 + np.sqrt(gamma_1d_2) * sig_2)), g_u]
H_int1_conj = [-0.5j * (a_u * (np.sqrt(gamma_1d_1) * sig_1.dag() + np.sqrt(gamma_1d_2) * sig_2.dag())), g_u.conjugate()]
gamma_phi_1 = 0.0 * gamma_1d_1
gamma_phi_2 = 0.0 * gamma_1d_2
Lo_us = [np.sqrt(gamma_1d_1)*sig_1, np.sqrt(gamma_1d_2)*sig_2, [a_u, g_u.conjugate()]]
# c_ops = [Lo_us, np.sqrt(gamma_phi_1) * sig_1.dag()*sig_1, np.sqrt(gamma_phi_2)*sig_2.dag()*sig_2]
c_ops = [Lo_us]
H = [H0, H_int1, H_int1_conj]
psi_0 = tensor(basis(emitter_dims, 0), basis(emitter_dims, 0), psi_coherent)
results = mesolve(H, psi_0, t, c_ops = c_ops)
ρ0 = psi_0 * psi_0.dag()
U = propagator(H, t, c_ops=c_ops) Code Output
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[10], line 2
1 ρ0 = psi_0 * psi_0.dag()
----> 2 U = propagator(H, t, c_ops=c_ops)
File /opt/anaconda3/envs/qutip5/lib/python3.11/site-packages/qutip/solver/propagator.py:66, in propagator(H, t, c_ops, args, options, **kwargs)
63 list_output = True
65 if not isinstance(H, (Qobj, QobjEvo)):
---> 66 H = QobjEvo(H, args=args, **kwargs)
68 if c_ops:
69 H = liouvillian(H, c_ops)
File /opt/anaconda3/envs/qutip5/lib/python3.11/site-packages/qutip/core/cy/qobjevo.pyx:223, in qutip.core.cy.qobjevo.QobjEvo.__init__()
File /opt/anaconda3/envs/qutip5/lib/python3.11/site-packages/qutip/core/cy/qobjevo.pyx:272, in qutip.core.cy.qobjevo.QobjEvo._read_element()
File /opt/anaconda3/envs/qutip5/lib/python3.11/site-packages/qutip/core/coefficient.py:175, in coefficient(base, tlist, args, args_ctypes, order, compile_opt, function_style, boundary_conditions, **kwargs)
173 for type_ in coefficient_builders:
174 if isinstance(base, type_):
--> 175 return coefficient_builders[type_](base, **kwargs)
177 if callable(base):
178 op = FunctionCoefficient(base, args.copy(), style=function_style)
File /opt/anaconda3/envs/qutip5/lib/python3.11/site-packages/qutip/core/cy/coefficient.pyx:418, in qutip.core.cy.coefficient.InterCoefficient.__init__()
ValueError: tlist must be the same len as the array to interpolate Expected Behaviour
Wanted it to return the propagator
Your EnvironmentQuTiP Version: 5.0.4 Numpy Version: 1.26.4 Scipy Version: 1.12.0 Cython Version: None Matplotlib Version: 3.9.2 Python Version: 3.11.9 Number of CPUs: 8 BLAS Info: Generic INTEL MKL Ext: False Platform Info: Darwin (x86_64)Additional Context
Even in Qutip 4.7, though the code runs, I'm seeing quite different outputs from mesolve and using the propagator (I might be using the propagator incorrectly though, so some help is appreciated):
plt.plot(t, expect(a_u.dag() * a_u, results.states), label="mesolve")
for i in range(len(t)):
plt.plot(t[i], expect(a_u.dag() * a_u, U[i]*ρ0), "x", color="orange")
plt.ylabel("Population")
plt.xlabel("Time")
plt.legend()
plt.show()
RetroSearch is an open source project built by @garambo | Open a GitHub Issue
Search and Browse the WWW like it's 1997 | Search results from DuckDuckGo
HTML:
3.2
| Encoding:
UTF-8
| Version:
0.7.4
