在Python QuTiP中集成不成功

I have been trying to use QuTiP to solve a quantum mechanics matrix differential equation (a Lindblad equation). Here is the code:

我一直在尝试用量子力学来解决量子力学矩阵微分方程(林德伯拉德方程)。这是代码:

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

This code is supposed to solve the following equation:

这段代码应该可以解决以下等式:

where L_i are matrices (in the list: [L1,L2,L3,L4,L5,L6,L7]), H is the hamiltonian, another matrix, is the density matrix, and is a constant equal to where T is the temperature, k is the Boltzmann constant, and , where h is Planck's constant.

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

Every time I run the code, it gives me the following error:

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

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.

After doing some debugging analysis, it seems like the first or second integration fails. The error tells me to increase the nsteps parameter, which I have tried. Even then it fails. Changing the list of times (the np.linspace function makes the list of times) also has no effect.

在做了一些调试分析之后，第一次或第二次集成似乎失败了。错误告诉我增加nsteps参数，我已经尝试过了。甚至失败。改变时间列表(np)linspace函数做时间列表)也没有效果。

I am desperate to know what I can do to fix this error. Please comment if you all need more details.

我非常想知道我能做些什么来修正这个错误。如果你们需要更多的细节，请评论。

Thanks for all your help!

谢谢你的帮助!

1 个解决方案

#1

From a programming point of view, the problem appears to be your value of gamma and therefore the size of your collapse operators. Print out gamma - it is of the order 10**25 - this seems to be what is preventing the solver from converging.

从编程的角度来看，问题似乎是您的伽玛值，因此是崩溃操作符的大小。打印输出伽玛-这是订单10* 25 -这似乎是阻止求解者收敛的原因。

Just for testing (I'm an engineer, not a quantum physicist...), I put in a smaller value of gamma (e.g. 0.1), the solver seems to work and gives apparently reasonable output in results.states

为了测试(我是工程师，不是量子物理学家……)，我输入了一个更小的伽玛值(例如0.1)，这个解算器似乎可以工作，并给出明显合理的结果

I don't quite understand your gamma - it seems to have units of cm^-1s^-2 as you have set it up. I wonder if you only want to divide by hbar once, maybe. As I say, I'm not a quantum physicist, so I'm only guessing here based on what makes the programming hang together combined with a bit of dimensional analysis.

我不太明白你的伽玛-它似乎有单位的cm-1 -2，你已经设置好了。我想知道你是否只需要除以hbar一次。就像我说的，我不是量子物理学家，所以我只是在猜测是什么让编程结合了一些维度分析。

EDIT

编辑

OP indicates in comments that the wrong order of magnitude / units for gamma does seem to be the programming issue (i.e. preventing numerical calculus from converging), but isn't totally clear on how to calculate gamma. At this stage, it may be worth posting a question at either http://physics.stackexchange.com or http://math.stackexchange.com about that - referencing this one for context if necessary.

OP在注释中指出，伽玛射线的大小/单位的错误顺序似乎确实是编程问题(即防止数值微积分收敛)，但对于如何计算伽玛射线并不完全清楚。在这个阶段，可能有必要在http://physics.stackexchange.com或http://math.stackexchange.com上发布一个问题——如果需要的话，可以在上下文中引用这个问题。

EDIT 2

编辑2

I note you asked this related question on the physics site. This makes it clear where the expression for gamma comes from and thereby clarifies that the constant terms presented as simply 30 and 150 in this question actually have units (Energy and frequency respectively). This changes the dimensional analysis - the units of gamma are s^-1 or, with appropriate conversion, cm^-1.

我注意到你在物理网站上问了这个相关的问题。这就清楚地说明了伽玛的表达式来自何处，从而澄清了在这个问题中仅仅表示为30和150的常数项实际上是有单位的(能量和频率各自)。这改变了维度分析——伽玛的单位是s-1，或者，通过适当的转换，是cm-1。

It also shows the value you mention in comments - 300 cm^-1.

它还显示了您在注释中提到的值——300 cm-1。