1

我一直在尝试使用 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 是哈密顿矩阵，另一个矩阵$\rho$是密度矩阵，并且$\伽马$是一个常数，$2\pi kT/\hbar*E_{R}/(\hbar\omega_{c})$其中 T 是温度, k 是玻尔兹曼常数, $\hbar$ = $h/2\pi$, 其中 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 函数生成时间列表）也没有效果。

我很想知道我能做些什么来解决这个错误。如果大家需要更多详细信息，请发表评论。

感谢你的帮助！

4

1 回答 1

1

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

只是为了测试（我是工程师，而不是量子物理学家......），我输入了一个较小的值gamma（例如 0.1），求解器似乎工作并给出了明显合理的输出results.states

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

编辑

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

编辑 2

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

它还显示了您在评论中提到的值 - 300 cm -1。

于 2015-06-03T01:00:23.283 回答