4. 仿真环境搭建与数值工具链
Lindblad主方程的数值仿真不是从写方程开始,而是从一条可信、可复现的工具链开始。本节要解决三个具体问题:如何把 Python 与 QuTiP 配置成能稳定运行的环境;如何正确调用 mesolve 与 mcsolve 这两个核心求解器;以及如何把仿真结果保存下来并做基础可视化。三者构成一条从"环境"到"计算"再到"数据"的完整链路,任何一环出错,后续的物理结论都不可信。
一、环境配置:版本锁定比安装本身更重要
QuTiP 的接口在不同大版本之间有实质性变化,尤其是 5.x 版本对数据后端和求解器返回值做了调整。因此配置的第一原则是锁定版本,而不是"装最新版"。推荐使用 conda 创建独立环境,把 Python、NumPy、SciPy、QuTiP 与绘图库一次性固定:
conda create -n lindblad python=3.11 -y
conda activate lindblad
conda install -c conda-forge qutip=5.0 numpy scipy matplotlib -y
python -c "import qutip; print(qutip.__version__)"
最后一行是必须执行的验证:能打印出版本号,说明 QuTiP 的 C 扩展编译正常。若在 Windows 上出现导入错误,多半是缺少与 Python 版本匹配的编译器或 conda 包通道混用,此时应优先用 conda-forge 单一通道重装,而不是逐个 pip 修补。环境验证还应包括一次最小物理计算,例如构造一个二能级系统的 sigmax() 并求本征值,确认线性代数后端可用。
二、求解器接口:mesolve 与 mcsolve 的分工
两个求解器面对的是同一物理模型,但回答的问题不同。mesolve 直接对密度矩阵做主方程演化,得到的是系综平均后的确定性格点;mcsolve 则通过量子轨迹抽样,用大量随机实现逼近同一平均结果,同时保留单次实现的随机性信息。选择依据是:需要平滑的平均值、期望值曲线时用 mesolve;需要研究单次测量记录、跳跃时刻分布或光子计数统计时用 mcsolve。
调用 mesolve 的关键参数是哈密顿量 H、初始态 rho0、时间数组 tlist、坍缩算符列表 c_ops,以及可选的期望值算符 e_ops。一个典型调用如下,其中 sigmam() 作为唯一耗散通道,模拟自发辐射:
import numpy as np, qutip as qt
H = 0.5 * 2*np.pi * qt.sigmax() # 拉比频率示例
psi0 = qt.basis(2, 0)
c_ops = [np.sqrt(0.1) * qt.sigmam()] # 衰减率示例
tlist = np.linspace(0, 20, 400)
res = qt.mesolve(H, psi0, tlist, c_ops,
e_ops=[qt.sigmax(), qt.sigmaz()])
print(res.expect[1][-1]) # 末时刻 <sigma_z>
这里 c_ops 中的系数是衰减率的平方根,这是 Lindblad 形式 L = sqrt(gamma) * sigma_minus 的直接体现,写错系数会系统性改变耗散强度。mcsolve 的接口几乎相同,但需要额外指定轨迹数 ntraj,且初始态通常用态矢量而非密度矩阵:
mc = qt.mcsolve(H, psi0, tlist, c_ops,
e_ops=[qt.sigmaz()], ntraj=500)
print(mc.expect[0][-1])
轨迹数 500 是示例值:轨迹太少,平均值噪声大;轨迹多则耗时线性增长。实践中可先用 100 条观察趋势,再用 1000 条以上做定量比较,并检查平均值是否随轨迹数增加而收敛到 mesolve 的结果——这是验证两条计算路径一致性的标准手段。
三、数据存储:让结果可复现、可比较
求解器返回的 Result 对象包含 states、expect 和 times 等字段,但直接 pickle 整个对象并不稳妥:QuTiP 对象结构随版本变化,跨环境读取容易失败。更可靠的做法是只保存数值数组与元数据。例如把时间轴、期望值、以及本次仿真的参数(哈密顿量系数、衰减率、轨迹数、随机种子)一起写入 .npz:
np.savez("run_gamma0.1.npz",
t=tlist,
sz=res.expect[1],
gamma=0.1, ntraj=0,
seed=1234, qutip=qt.__version__)
把参数与版本号写进文件,是为了让半年后的自己或合作者能判断这条曲线是在什么条件下产生的。若需要人工查看,可另存一份 CSV;若数据量大,则按"一次扫描一个文件"的方式命名,避免单文件膨胀。随机种子对 mcsolve 尤其重要,固定种子才能复现同一组轨迹。
四、可视化基础:先看对,再看美
基础可视化只需回答两个问题:期望值是否收敛到合理终值,以及耗散是否按预期速率衰减。对上面的二能级例子,⟨σz⟩ 应从初值 +1 单调趋向 -1 附近,衰减时间尺度约由 γ 决定。绘图时把 mesolve 与 mcsolve 的平均曲线画在同一坐标系,若两者在轨迹数足够时重合,说明模型与求解设置一致;若明显偏离,应优先检查 c_ops 系数、初态类型(态矢量与密度矩阵混用)以及 tlist 是否覆盖了完整演化过程。
check_env.py 脚本,每次开始新仿真前运行:打印 QuTiP 版本、跑一次 10 步的 mesolve 小例子、确认 np.savez 可写。三秒钟的检查能避免大量"结果异常其实是环境坏了"的排查时间。