python - Python QuTiP 中的集成未成功

Question

我一直在尝试使用 QuTiP 来求解量子力学矩阵微分方程(Lindblad 方程)。这是代码:

from qutip import *
from matplotlib import *
import numpy as np

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)
ground=Qobj(np.array([[ 0.0863685 ],
  [ 0.17141713],
  [-0.91780802],
  [-0.33999268],
  [-0.04835763],
  [-0.01859027],
  [-0.05006013]]))

rho0 = ground*ground.dag()
from scipy.constants import *
ktuple=physical_constants['Boltzmann constant in eV/K']
k = ktuple[0]* 8065.6
htuple = physical_constants['Planck constant in eV s'] 
hbar = (htuple[0]* 8065.6)/(2*pi)
gamma=(2*pi)*((k*300)/hbar)*(35/(150*hbar))

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=100000)

results = mesolve(H, rho0, np.linspace(0.0, 1000, 20),[L1,L2,L3,L4,L5,L6,L7],[],options=options)
print results

这段代码应该求解以下方程:

其中 L_i 是矩阵(在列表中:[L1,L2,L3,L4,L5,L6,L7])，H 是哈密顿量，另一个矩阵，是密度矩阵，是一个常数，等于 其中 T 是温度，k 是玻尔兹曼常数， ，其中 h 是普朗克常数。

每次运行代码时，都会出现以下错误:

ZVODE--  At T (=R1) and step size H (=R2), the
       corrector convergence failed repeatedly
       or with abs(H) = HMIN
      In above,  R1 =  0.0000000000000D+00   R2 =  0.1202322246215D-36
/usr/local/lib/python2.7/dist-packages/scipy/integrate/_ode.py:853: UserWarning: zvode: Repeated convergence failures. (Perhaps bad Jacobian supplied or wrong choice of MF o
r tolerances.)
  'Unexpected istate=%s' % istate))
Traceback (most recent call last):
  File "lindbladqutip.py", line 48, in <module>
    results = mesolve(H, rho0, np.linspace(0.0, 1000, 20),[L1,L2,L3,L4,L5,L6,L7],[],options=options)
  File "/projects/d6138712-e5f4-4d85-9d4d-77ce0a7b4a61/.local/lib/python2.7/site-packages/qutip/mesolve.py", line 264, in mesolve
    progress_bar)
  File "/projects/d6138712-e5f4-4d85-9d4d-77ce0a7b4a61/.local/lib/python2.7/site-packages/qutip/mesolve.py", line 692, in _mesolve_const
    return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)
  File "/projects/d6138712-e5f4-4d85-9d4d-77ce0a7b4a61/.local/lib/python2.7/site-packages/qutip/mesolve.py", line 866, in _generic_ode_solve
    raise Exception("ODE integration error: Try to increase "
Exception: ODE integration error: Try to increase the allowed number of substeps by increasing the nsteps parameter in the Options class.

经过一些调试分析，似乎第一次或第二次集成失败了。错误告诉我增加nsteps参数，我已经尝试过了。即使如此，它还是失败了。更改时间列表(np.linspace 函数生成时间列表)也没有效果。

我非常想知道我能做些什么来修复这个错误。如果您需要更多详细信息，请发表评论。

感谢您的帮助!

Answer 1

从编程的角度来看，问题似乎出在 gamma 的值上。以及折叠操作符的大小。打印出gamma - 它的顺序是10**25 - 这似乎是阻止求解器收敛的原因。

只是为了测试(我是工程师，不是量子物理学家......)，我输入了较小的值 gamma (例如 0.1)，求解器似乎可以工作，并在 results.states 中给出明显合理的输出

我不太明白你的gamma - 正如您所设置的，它似乎的单位为cm^-1s^-2。我想知道你是否只想除以 hbar一次，也许。正如我所说，我不是量子物理学家，所以我只是根据使编程结合在一起的原因并结合一些维度分析来猜测。

编辑

OP 在评论中指出 gamma 的数量级/单位错误。似乎确实是编程问题(即防止数值演算收敛)，但并不完全清楚如何计算 Gamma 。在此阶段，可能值得在 http://physics.stackexchange.com 中发布问题。或http://math.stackexchange.com关于这一点 - 如有必要，请引用此内容以获取上下文。

编辑 2

我注意到你问了这个related question on the physics site 。这清楚地表明 the expression for gamma comes from 在哪里从而澄清了常数项简单地表示为 30和150在这个问题中实际上有单位(分别是能量和频率)。这改变了维度分析 - Gamma 的单位为 s^-1，或者通过适当的转换为 cm^-1。

它还显示您在评论中提到的值 - 300 厘米^-1。