我已经为 open quantum system 建模工作了很长时间使用 Lindblad Equation .哈密顿量如下:
但是,哈密顿量中添加了另外两个矩阵。其中之一的所有对角线项都等于 -33.3333i,其他所有项都为零。另一个是第三个对角线项等于 -0.033333i 的矩阵。
Lindblad 方程是这样的:
其中 L_i 是矩阵(在列表中:[L1,L2,L3,L4,L5,L6,L7])。 L_i 的矩阵只是一个 7x7 矩阵,除 L_(ii)=1 外全为零。 H 是总哈密顿量,是密度矩阵,是等于 的常数其中 T 是温度,k 是玻尔兹曼常数, ,其中 h 是普朗克常数。 (注意 gamma 在 natural units 中)
以下代码求解 Lindblad 方程,从而计算密度矩阵。然后计算并绘制此与时间的关系:
这被称为站点 3 人口。 被称为胸罩和被称为凯特。两者都是载体。在这种情况下,请参阅代码以了解它们的定义。
代码如下:
from qutip import Qobj, Options, mesolve
import numpy as np
import scipy
from math import *
import matplotlib.pyplot as plt
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]
])
recomb = np.zeros((7, 7), dtype=complex)
np.fill_diagonal(recomb, 33.33333333)
recomb = recomb * -1j
trap = np.zeros((7, 7), complex)
trap[2][2] = -0.033333333333j
hamiltonian = recomb + trap + hamiltonian
H = Qobj(hamiltonian)
# Note the extra .0 on the end to convert to float
gamma = (2 * pi) * (296 * 0.695) * (35.0 / 150)
L1 = 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 = 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 = 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 = 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 = 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 = 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 = 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:
rho0=Qobj(L1)
L1 = Qobj(gamma * L1)
L2 = Qobj(gamma * L2)
L3 = Qobj(gamma * L3)
L4 = Qobj(gamma * L4)
L5 = Qobj(gamma * L5)
L6 = Qobj(gamma * L6)
L7 = Qobj(gamma * L7)
options = Options(nsteps=1000000, atol=1e-5)
bra3 = [[0, 0, 1, 0, 0, 0, 0]]
bra3q = Qobj(bra3)
ket3 = [[0], [0], [1], [0], [0], [0], [0]]
ket3q = Qobj(ket3)
starttime = 0
# this is effectively just a label - `mesolve` alwasys starts from `rho0` -
# it's just saying what we're going to call the time at t0
endtime = 100
# Arbitrary - this solves with the options above
# (max 1 million iterations to converge - tolerance 1e-10)
num_intermediate_state = 100
state_evaluation_times = np.linspace(
starttime,
endtime,
num_intermediate_state
)
result = mesolve(
H,
rho0,
state_evaluation_times,
[L1, L2, L3, L4, L5, L6, L7],
[],
options=options
)
number_of_interest = bra3q * (result.states * ket3q)
points_to_plot = []
for number in number_of_interest:
if number == number_of_interest[0]:
points_to_plot.append(0)
else:
points_to_plot.append(number.data.data.real[0])
plt.plot(state_evaluation_times, points_to_plot)
plt.show()
exit()
此代码使用名为 qutip 的 Python 模块.它有一个使用 scipy.integrate.odeint 的内置 Lindblad 方程求解器。
目前,这个程序显示这个:
然而,站点 3 的人口上限应该为 0。因此,它应该会缓慢下降到零。特别是到 t=75 时,应该开始减少。
这段代码运行了,但没有像我解释的那样产生正确的结果。那么现在,为什么它没有产生正确的结果呢?我的代码有问题吗?
我查看了我的代码,查看每一行是否与我正在使用的模型相匹配。他们完美匹配。问题一定出在代码上,而不是物理上。
我做了一些调试提示,所有的矩阵和 Gamma 都是正确的。但是,我仍然怀疑 trap
矩阵中有什么东西。我这么认为的原因是因为情节看起来像系统的动力学没有trap
矩阵,我的陷阱矩阵的定义会不会有问题没注意到?
请注意,代码需要几分钟才能运行。运行代码时请耐心等待!
最佳答案
(注意:这是我希望在编程意义上的答案,而不是物理答案。)
我独立运行了你的模拟,没有使用 qutip
,我得到的结果基本相同。所以好消息(也许?:))是这不是您的编程问题,而是物理问题或至少是“选择参数”问题。
这是我的结果:
还有一个我工作的笔记本,参数都和你的一样,除了不同的时间尺度(下面解释)。我使用与 qutip
相同的集成方法但不是 qutip 本身:Notebook Link .
一些注意事项:
当您执行
from math import *
你导入一个函数gamma
, 然后你命名一个变量gamma
,这给我带来了麻烦,你以后可能要小心了。当您将 linblad 运算符乘以
\gamma
时而不是总和,它们将在主方程中出现两次,因此您实际上是在指定\gamma^2
这里。这会影响时间尺度。<3|rho(t)|3>
只是第三个对角矩阵元素,你在这里真的不需要内积。
物理/参数方面需要检查的几件事。
- 从您链接的论文中,
- 是
\Gamma
肯定是 100/3? - 是
\kappa_3
肯定是 0.1/3 而所有其他人是 0? - 初始状态肯定是所有种群都处于0状态吗?
- 是
- 我不是最新的能量转移模型,但这里的哈密尔顿是非厄尔米特的,非平凡的虚部(尽管仍然很小)是在密度矩阵对角线上生成的。确保您确切了解这些人如何以及为何使用此模型,因为我觉得它很奇怪!
关于python - 开放量子系统建模,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/34669082/