4. 量子系统仿真环境搭建

参数扫描与敏感性分析的前提，是有一个"参数可调、结果可复现"的仿真环境。如果每次改参数都要手工改代码、跑一次、抄一次结果，那么扫描几十个参数点就会变成体力活，而且极易出错。本节的目标是把量子系统封装成参数化模型，再用脚本把"参数网格 → 批量求解 → 结果汇总"这条链路自动化。工具上以 QuTiP 为主，Qiskit 作为补充说明。

一、环境准备与版本约束

QuTiP 依赖 NumPy、SciPy、Cython，建议用独立虚拟环境安装，避免与 Qiskit 的依赖冲突。一个可用的最小环境如下：

python -m venv qenv
source qenv/bin/activate        # Windows: qenv\Scripts\activate
pip install "qutip>=5.0" numpy scipy matplotlib
pip install qiskit qiskit-aer   # 需要对比或做门级仿真时安装
python -c "import qutip; print(qutip.__version__)"

QuTiP 5 起底层数据结构有较大变化，Qobj 的稀疏存储与求解器接口与 4.x 不完全一致。课程示例统一按 QuTiP 5 的写法给出；若你的环境是 4.x，mesolve 的调用方式基本兼容，但要注意 qutip.solver 命名空间的变化。版本确认这一步不要跳过，很多"结果对不上"的问题其实来自版本差异。

二、把系统写成参数化函数

参数化建模的核心思想是：系统的一切可变部分都作为函数入参，函数内部只负责"由参数构造哈密顿量与耗散算符"，不负责扫描和保存。以驱动二能级系统（Rabi 模型）为例，设驱动强度为 Ω、失谐为 Δ、弛豫率 γ、退相位率 γφ，则：

import numpy as np
import qutip as qt

sx, sy, sz = qt.sigmax(), qt.sigmay(), qt.sigmaz()
sm = qt.destroy(2)          # 降算符 |0><1|

def build_system(Omega, Delta, gamma, gamma_phi):
    """返回 (H, c_ops)，所有物理量均为角频率单位"""
    H = 0.5 * Omega * sx + 0.5 * Delta * sz
    c_ops = [np.sqrt(gamma) * sm,
             np.sqrt(gamma_phi) * sz]
    return H, c_ops

def run_point(Omega, Delta, gamma, gamma_phi, tlist, psi0):
    H, c_ops = build_system(Omega, Delta, gamma, gamma_phi)
    return qt.mesolve(H, psi0, tlist, c_ops,
                      e_ops=[sz, sm.dag() * sm])

这里有两个约定必须先说清楚，否则后续所有图和数值都会失去意义。第一，单位约定：Ω、Δ、γ 全部用角频率（rad/s 或无量纲的 1/time），时间轴 tlist 用同一时间单位，这样 Omega * t 才是无量纲相位。第二，算符约定：sz 的本征值为 ±1，sm.dag()*sm 对应激发态布居，取值 0～1。参照约定一旦定下，全文与图示保持一致。

把 build_system 与 run_point 分开，是为了让"物理模型"和"求解流程"解耦：换一个模型（如三能级、腔 QED）只需重写 build_system，扫描脚本完全不用动。

三、单参数扫描：从网格到曲线

单参数扫描最典型的任务是"固定其他参数，看某一观测量随一个参数如何变化"。例如固定 Δ=0、γ=0.05、γφ=0.02，让 Ω 从 0.1 扫到 2.0，共 20 个点，观察 t=10 时刻的激发态布居：

tlist = np.linspace(0, 10, 201)
psi0 = qt.basis(2, 0)                 # 初态 |0>
Omegas = np.linspace(0.1, 2.0, 20)
pe = []
for Om in Omegas:
    r = run_point(Om, 0.0, 0.05, 0.02, tlist, psi0)
    pe.append(r.expect[1][-1])        # 末时刻激发态布居
pe = np.array(pe)

结果怎么读？在弱驱动、短时间近似下，激发态布居约正比于 (Ω/√(Ω²+Δ²))²·sin²(Ω_eff·t/2)，因此 Ω 增大时末时刻布居会呈现非单调的振荡式上升——这是"固定观测时刻"带来的相位效应，不是物理上的单调响应。判断扫描是否合理，可以先看曲线是否光滑、是否出现明显跳变；若相邻点差异极大，往往是 tlist 分辨率不足或观测时刻恰好落在振荡节点上，应加密时间网格或改用时间平均量（如 np.mean(r.expect[1])）作为扫描指标。

四、多参数扫描：网格、并行与结果组织

多参数扫描的本质是笛卡尔积。以 Ω 与 Δ 两个参数为例，用 itertools.product 生成网格，并把每个点的结果连同参数一起写入结构化记录，避免"结果和参数对不上"：

import itertools, json

records = []
for Om, De in itertools.product(np.linspace(0.1, 2.0, 10),
                                np.linspace(-1.0, 1.0, 9)):
    r = run_point(Om, De, 0.05, 0.02, tlist, psi0)
    records.append({"Omega": float(Om), "Delta": float(De),
                    "pe_final": float(r.expect[1][-1]),
                    "pe_mean": float(np.mean(r.expect[1]))})
with open("scan_2d.json", "w") as f:
    json.dump(records, f, ensure_ascii=False, indent=1)

当网格点数达到数百以上时，串行循环会明显变慢。QuTiP 的 mesolve 本身是单线程的，可行的加速方式是用 multiprocessing.Pool 按参数点并行，每个进程独立构造 Qobj（注意不要在进程间共享 Qobj，稀疏矩阵序列化开销大）。并行时务必固定随机种子（若模型含随机项），并保证每个任务只依赖自己的参数，这样结果才可复现。

结果组织上，二维扫描建议直接整理成矩阵再画热图：行对应 Ω、列对应 Δ，单元格填 pe_final。热图上若出现一条沿 Δ=0 的亮带，说明共振附近响应最强；若亮带随 Ω 增大而变宽，则提示驱动展宽效应。这类"先看整体形状、再回到单点核对"的读图顺序，比逐点看数值高效得多。

五、Qiskit 路线的差异

若系统本身是门级线路模型，Qiskit 更合适：用 QuantumCircuit 构造参数化线路，用 Parameters 绑定数值，再交给 AerSimulator 或 transpile 后的后端执行。与 QuTiP 的关键差异在于：Qiskit 扫描的是线路参数与采样次数，输出是计数分布而非连续期望值，因此敏感性指标通常取某个比特的 P(1) 或期望值 ⟨Z⟩，并需要足够多的 shots 来压低统计噪声。做敏感性分析时，shots 不足会让"参数微小变化引起的响应差异"淹没在采样涨落里，这一点在 QuTiP 的确定性求解中不存在。

六、一个可直接复用的扫描配置模板

把扫描任务写成配置，而不是散落在脚本里的字面量，是让实验可复现的关键。下面这份配置表可直接作为课程练习的起点，字段只保留与本节任务真正相关的项：

至此，环境搭建的交付物应当是三样东西：一个只依赖参数的 build_system、一个统一的 run_point 求解入口、一份把参数网格与结果记录绑定在一起的扫描配置。有了这三样，后续的单参数曲线、二维热图乃至敏感性指标计算，都只是在同一套流水线上更换观测量与参数范围而已。