量子相位估计QPE这个词第一次看到它出现在分子动力学语境里我第一反应是这俩不是一回事。量子化学求的是定态能量分子动力学跑的是原子核轨迹中间隔着玻恩-奥本海默近似和一整条能量梯度链。但真把流程捋顺之后会发现QPE 恰好卡在整条链路上最关键、也最贵的那个环节单点能量的高精度求解。用 Qiskit 把 QPE 搭出来再接到经典积分器上就能拼出一条量子算能量、经典算轨迹的混合路线。这篇东西写给两类人一类是已经会跑简单的量子线路模拟、想往化学方向走的人一类是做经典分子动力学做烦了、想看看量子这条腿到底能补多少精度的人。下面我把从哈密顿量构造到电路实现、从参数推导到踩坑记录的全过程摊开讲代码都是能直接粘进 notebook 跑的。1. 先把边界划清QPE 在分子动力学里算的到底是哪一段1.1 经典 MD 的精度天花板卡在势能面这一步做经典分子动力学的人都清楚轨迹本身不难跑速度 Verlet、蛙跳积分、周期性边界、恒温器这些东西几十年了没什么大变化。真正决定你结果靠不靠谱的是你每一步问的那个力场值准不准。力场这条路走的是预先拟合把大量高精度量化计算的结果拟合成功率场参数之后每次调用只做多项式求值。代价是换体系就得重新拟合而且对那些成键断裂、过渡态、强关联电子体系拟合出来的势能面根本没法用。另一条路叫从头算分子动力学跑一步算一次电子结构通常用密度泛函理论。这条路精度确实上去了但成本是随电子数三到四次方甚至更高往上翻。一个几十个原子的体系跑几百飞秒的轨迹能把一个集群排满算上几个月。所以这些年的实际做法全是在中间找折中小体系用高精度方法大体系用泛函加修正。量子计算介入的价值就在这里——它提供的是一条理论上可扩展、精度可控的路径不是靠拟合、不是靠近似交换关联泛函而是直接把哈密顿量的本征值读出来。1.2 QPE 输出的是本征能量不是轨迹这一点必须掰开说清楚因为很多人第一眼会误解。QPE 给定一个西算子 U 和它的一个本征态 |ψ⟩返回的是 U 的本征相位 φ满足U|ψ⟩ e^(2πiφ)|ψ⟩在量子模拟场景里我们把分子哈密顿量 H 编码成时间演化算子 U e^(iHt)那么U|ψ⟩ e^(iEt)|ψ⟩相位和能量之间就只差一个时间因子φ Et/(2π)。所以 QPE 干的事情是把一个能量值翻译成一个相位值再用干涉的办法把这个相位读成二进制小数。它给不了你原子坐标怎么变给不了你温度、压力、扩散系数。它给的是一个数在当前核构型下这个电子体系的能量是多少。要做动力学还得你自己在外面套一层经典积分器对每个核构型调一次 QPE 拿单点能再数值求梯度算力再积分。这就是我说的混合流程。1.3 我采用的混合流程量子内核负责单点能具体跑起来大概是这样一条链用经典方法分子力场或者低精度泛函先跑一段预平衡拿到一批核构型的采样点对每个关键构型构造对应的电子哈密顿量映射到量子比特上 QPE 求基态能量精度目标是化学精度也就是大约 1.6 mHa1 kcal/mol 量级把这些高精度单点能作为锚点反过来校正或拟合低成本方法得到一条修正后的势能面用修正后的势能面继续跑 MD。注意第 4 步我这里没有天真到每个 MD 步都调一次 QPE——那在当前硬件上完全不现实。真正能落地的用法是稀疏锚点校正QPE 负责提供少量高置信度的参考值把经典低精度方法的系统偏差压下去。这个思路和目前在用的 Δ-machine learning 是一脉相承的只不过参考值的来源换成了量子计算。提示如果你看到有人写用量子计算加速分子动力学先问一句他的 QPE 是每步调用还是抽样调用。前者目前基本是纸上方案后者才有工程可行性。1.4 工具链为什么落在 Qiskit 上选型这件事我也纠结过。Qiskit 的优势不在单个模块的性能而在于它把三样东西凑齐了电路构建层QuantumCircuit、各种 standard gate、化学模块qiskit-nature 提供从原子坐标到二次量子化哈密顿量的完整链路、以及本地模拟器Aer 的 statevector 和 matrix product state 后端。对比一下常见的几种组合有的框架在化学映射这块更成熟但电路层的可控性差一些有的模拟器性能强但从费米子到比特的映射得自己写。Qiskit 在这三块的水平都在够用且顺手的区间尤其是PauliEvolutionGate和它的.control()方法省掉了我手写受控 Trotter 电路的一大堆麻烦。另外SparsePauliOp和numpy之间的来回转换很干净做完 QPE 想拿经典对角化结果做交叉验证一行代码的事。至于后端本地就用 Aer 的 statevector 模拟器起步需要跑更大比特数再考虑 MPS 方法或者真实硬件。这里有个取舍要提前想清楚statevector 的内存是 2^n 复振幅n20 大概就 16 GB 起步了n25 基本没戏。所以本地验证阶段把比特数控制在 20 以内超出的部分走 MPS 或者分布式。2. 开跑之前环境、比特数与演化时间三件套2.1 依赖版本与最小可跑环境先说我用的版本组合这套是跑通的。Qiskit 从 1.0 开始做了不小的 API 调整很多老教程里的qiskit.opflow已经废弃了直接照抄会报一堆 deprecation warning 甚至 ImportError所以版本必须锁。pip install qiskit1.2,1.3 qiskit-aer0.15 qiskit-nature0.7,0.8 pyscf三块东西各管一摊qiskit负责电路和算符表达qiskit-aer提供模拟后端qiskit-nature处理化学部分从原子坐标生成哈密顿量。pyscf是底层做积分计算的qiskit-nature的驱动会调它。import 部分我习惯这样写能少踩很多坑import numpy as np from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile from qiskit.circuit.library import QFT, PauliEvolutionGate from qiskit.quantum_info import SparsePauliOp from qiskit.synthesis import LieTrotter from qiskit_aer import AerSimulator注意PauliEvolutionGate和LieTrotter这两个类的位置在不同小版本间挪过如果 import 失败先去查你装的那个版本的qiskit.synthesis下有没有LieTrotter。别急着怀疑自己写错了。2.2 相位比特数怎么倒推出来这是整个方案里第一个要算清楚的东西因为它直接决定你的电路有多大、能跑到什么精度。记 QPE 的相位寄存器有 n 个比特总演化时间最大幂次那一路是 t。相位被表示成 n 位二进制小数最小可分辨相位是 1/2^n。换算成能量ΔE 2π / (t · 2^n)这个式子是所有参数选择的出发点。反过来说给定你需要的能量精度 ΔE 和选定的演化时间 t比特数就是n ceil( log2( 2π / (t · ΔE) ) )举个有实际意义的算例。假设我们的目标体系是氢分子用 STO-3G 基组经过活性空间约化后本征值跨度大约 2.5 Ha从基态到较高激发态。如果取 t 2π/2.5 ≈ 2.51那么n 6ΔE 2.5 / 64 ≈ 39 mHa只能看个大概n 10ΔE 2.5 / 1024 ≈ 2.44 mHa接近化学精度边缘n 11ΔE ≈ 1.22 mHa达到 1 kcal/mol 以内的要求n 14ΔE ≈ 0.15 mHa可以讨论色散力级别的差异了。所以在本地模拟能承受的比特预算下加上目标寄存器的比特数总电路宽度最好控制在 20 以内n 取 10 到 12 是个甜点区。对应电路宽度大概是 12 4 16 比特statevector 模拟毫无压力跑几千 shots 也就几秒钟。2.3 演化时间 t 的取值推导附算例t 不是一个随便填的数它被两头夹着。下界来自精度需求t ≥ 2π / (ΔE · 2^n)你把上一节的式子反过来看就明白了t 越大同样的比特数能分辨的能量间隔越细。上界来自能量窗口QPE 读出的 φ 是模 1 的也就是只能区分相差 2π/t 整数倍的能量。要想让解唯一必须保证整个能量跨度小于这个周期t 2π / (E_max - E_min)这两个约束一夹t 的最优选择其实很明确取窗口宽度的倒数乘上 2π再留 10% 到 20% 的余量。t 2π / (1.2 · (E_max - E_min))用上一节的 H2 算例跨度 2.5 Hat 2π / (1.2 × 2.5) 6.2832 / 3.0 ≈ 2.094这样能量窗口实际覆盖 2π/2.094 ≈ 3.0 Ha比真实跨度 2.5 Ha 宽出 0.5 Ha 的缓冲既不会因为估错跨度导致卷绕也不会浪费精度。提示不要为了追求精度盲目把 t 调大。t 翻倍能量窗口就减半一旦真实能量超出窗口你得到的结果是经过卷绕的假值而且它看起来完全正常——这是最难排查的一类错误。2.4 能量卷绕最容易翻车的地方接着上一节说。QPE 输出 φ ∈ [0, 1)对应的能量是E 2πφ/t 2πk/t k 为整数这个2πk/t就是别名一串等间隔的候选能量你没法从单次测量里区分它们。解决思路有三条第一条是平移哈密顿量。把 H 整体减去一个保守下界 E_ref让所有本征值落在 [0, 2π/t) 区间内。E_ref 可以取经典方法算出来的基态能量再往下压个 1~2 Ha保守一点没关系只是稍微浪费一点窗口宽度。平移之后 k 0直接读 φ 就行。第二条是用先验知识掐掉多余的解。你知道基态能量的大致范围比如经典 DFT 给了个 ±0.5 Ha 的估计那就只接受落在这个范围内的 k。第三条是多组 t 值交叉验证。用两个互质的演化时间各跑一次只有真实能量会同时满足两组方程。这招在实验物理里是老套路了代价是电路要跑两遍。我在代码里默认用第一条简单可靠。具体写法是把哈密顿量加一个常数项E_ref -2.0 # 保守下界一定要比真实基态低 identity_shift SparsePauliOp.from_list([(I * num_qubits, -E_ref)]) H_shifted H identity_shift注意I * num_qubits这个写法SparsePauliOp 的 Pauli 字符串长度必须和比特数严格一致少一个字符就会抛异常而且报错信息不太直观。3. 分子哈密顿量落地的全过程3.1 从原子坐标到费米子算符这一段是整个流程里最容易被人跳过、但出错代价最高的环节。你手上有的是一串原子坐标和原子序数QPE 需要的是一个作用在量子比特上的厄米算符。中间隔着两层转换。第一层是玻恩-奥本海默近似下的电子哈密顿量H -Σ_i ∇_i²/2 - Σ_i Σ_A Z_A/|r_i - R_A| Σ_{ij} 1/|r_i - r_j| Σ_{AB} Z_A Z_B/|R_A - R_B|四项分别是电子动能、电子-核吸引、电子-电子排斥、核-核排斥。最后一项是个常数在固定核构型下可以先算出来存着最后再加回去。第二层是二次量子化。把上面这个算符用一组单电子基函数通常就是原子轨道线性组合成的分子轨道展开得到H Σ_pq h_pq a_p† a_q (1/2) Σ_pqrs h_pqrs a_p† a_q† a_r a_s其中h_pq是单电子积分h_pqrs是双电子积分。这些积分本身是经典可算的qiskit-nature底下的驱动会调pyscf完成这一步。整个过程对使用者是透明的但你要知道它在算什么否则后面比特数对不上、能量差了核排斥项这些问题都会让你抓瞎。3.2 Jordan-Wigner 和 Bravyi-Kitaev 该怎么挑费米子算符满足反对易关系量子比特上的 Pauli 算符不满足所以中间必须插一层映射。两种主流选择Jordan-Wigner 映射最直观一个费米子模式对应一个比特占据数为 1 对应 |1⟩。代价是引入一串 Z 算符组成的尾巴来维持反对易性导致 Pauli 串的长度随比特序号线性增长。对 4 到 8 比特的小体系问题不大比特数一多长 Pauli 串会让 Trotter 电路的深度爆炸。Bravyi-Kitaev 映射用一套更巧妙的编码把占据信息分散到多个比特上Pauli 串长度降到对数级。代价是线路结构不直观调试的时候很难从比特状态直接读出哪个轨道被占了。另外它要求比特数是 2 的幂次超出的部分要靠补充虚拟比特处理。我的选择是比特数 ≤ 10 用 JW超过 10 用 BK 或者 Parity。做教学演示和结果验证的时候 JW 的直观性价值很高你能一眼看出某个比特翻转对应哪个自旋轨道大规模跑的时候 BK 的线路深度优势就压过直观性了。还有第三个选项 Parity 映射它把粒子数和自旋对称性直接编码进去配合两比特约化能把 H2 这类小分子压到 2 个比特。这个后面单独讲。3.3 活性空间约化先砍比特再谈精度一个水分子用 STO-3G 基组是 7 个原子、14 个电子、14 个空间轨道映射到量子比特是 28 个。这个宽度在当前模拟器上勉强能跑 statevector但受控时间演化那部分会生成海量的 CNOT模拟时间直接上天。活路是活性空间近似把能量最低的芯电子冻结掉它们基本不参与化学反应把能量最高的虚轨道扔掉它们贡献很小只保留中间的活性轨道。一个常用的经验是保留最高占据轨道能量 - 10 eV到最低空轨道能量 10 eV这个窗口内的轨道。对 H2 分子2 个电子、2 个空间轨道JW 映射后是 4 个比特。如果再叠加 Parity 映射的对称性约化能压到 2 个比特——这就是为什么所有 QPE 教程都爱拿 H2 举例它小到可以在笔记本上跑完整个流程。但别被这个骗了。活性空间是近似砍掉的部分不参与关联能计算所以你算出来的能量会系统性地高于真实值。做方法验证没问题做实际精度评估的时候要把这部分误差算进预算里。3.4 实操搭出 qubit 哈密顿量并做经典校验真实分子的路径是这样from qiskit_nature.second_q.drivers import PySCFDriver from qiskit_nature.second_q.mappers import JordanWignerMapper driver PySCFDriver(atomH 0 0 0; H 0 0 0.735, basissto3g) problem driver.run() mapper JordanWignerMapper() fermionic_op problem.hamiltonian.second_q_op() qubit_op mapper.map(fermionic_op) # 别忘了核排斥能它是个经典常数 nuclear_repulsion problem.hamiltonian.nuclear_repulsion_energy print(qubit_op)nuclear_repulsion这一项单独拎出来是有原因的它不进量子线路就是个常数最后加回到能量结果上。我第一次跑的时候忘了加算出来的能量比参考值差了 0.7 Ha排查了整整一个下午。为了把整条流程验证通我建议先用一个手写的小哈密顿量做通关可控性高出问题好定位num_qubits 2 H SparsePauliOp.from_list([ (II, 0.0), (ZI, -0.4), (IZ, -0.4), (XX, 0.15), (YY, 0.15), (ZZ, 0.1), ]) # 经典对角化拿参考值 eigs np.linalg.eigvalsh(H.to_matrix()) print(精确本征值:, np.round(eigs, 6))这两个比特的哈密顿量在结构上模拟了一个简化版的电子关联体系两个局域能量项、一个 ZZ 相互作用、加上 XX 和 YY 组成的交换/配对项。它的基态是四个计算基态的叠加正好可以考察 QPE 在初态不是本征态时的行为。注意SparsePauliOp.from_list里的 Pauli 字符串比特顺序是左高右低也就是第一个字符对应编号最大的比特。这点和某些其他工具的约定相反跨工具对比结果的时候务必先确认否则你会得到一组看起来对得上、实际下标全反的结果。4. QPE 电路的手工搭法4.1 电路骨架拆成四段看标准 QPE 的线路分四段我按执行顺序捋第一段相位寄存器初始化。n 个相位比特全部置 |0⟩然后每个加一个 Hadamard。这一步的作用是制造叠加态让这 n 个比特同时探测所有可能的相位分量。第二段受控时间演化。对第 j 个相位比特j 从 0 到 n-1施加受控的 U^(2^j)。也就是 U 作用 1 次、2 次、4 次……一直到 2^(n-1) 次。这是把相位信息踢到相位寄存器上的关键步骤等效于做了一次相位反冲phase kickback。第三段逆量子傅里叶变换。第二段结束后相位寄存器处于一个傅里叶基下的态能量信息藏在相对相位里。逆 QFT 把它转回计算基使测量结果直接对应相位的二进制表示。第四段测量。在计算基上测相位寄存器得到 n 位比特串。代码骨架大致是qr_e QuantumRegister(n_eval, eval) qr_t QuantumRegister(num_qubits, tgt) cr ClassicalRegister(n_eval, c) qc QuantumCircuit(qr_e, qr_t, cr) qc.h(qr_e) # 第一段 for j in range(n_eval): # 第二段 ... qc.append(QFT(n_eval, inverseTrue, do_swapsTrue).to_gate(), qr_e) # 第三段 qc.measure(qr_e, cr) # 第四段4.2 受控时间演化算子的构造细节这是整个电路里最重的一块。U^(2^j) 不能真的把一个酉矩阵自乘 2^j 次再求受控版本——那矩阵维度直接爆。正确做法是把这个幂次吸收到时间参数里U^(2^j) e^(iH·t·2^j)也就是对第 j 个比特用演化时间t·2^j构造同一个哈密顿量的演化门然后取受控。Qiskit 里PauliEvolutionGate会自动把哈密顿量分解成 Pauli 串的和每一项单独做指数化再按 Trotter 公式拼接。用法for j in range(n_eval): evo PauliEvolutionGate( H_shifted, time-t * (2 ** j), # 负号实现 e^{iHt} synthesisLieTrotter(reps2) # 二阶 Trotter 分解 ) ctrl_evo evo.control(1) # 加一个控制比特 qc.append(ctrl_evo, [qr_e[j]] list(qr_t))关于time-t这个负号值得单拎出来说。PauliEvolutionGate默认实现的是e^(-iHt)而我们希望 U 的本征相位满足φ Et/(2π)直接对应正的能量值。取负时间等价于实现e^(iHt)省去了事后在换算公式里加符号的麻烦。你也可以不用负号靠公式E -2πφ/t反解两者完全等价只是前者读代码时不容易搞混。LieTrotter(reps2)里的 reps 是 Trotter 步数。一阶 Trotter 的误差是 O(t²/reps)二阶是 O(t³/reps²)。对小体系、小演化时间reps1 到 2 足够了。判断标准是把 reps 从 2 加到 4如果结果变化小于你要求的精度就不用再加了。提示.control(1)返回的门的比特顺序是控制比特在前目标比特随后。append 的时候如果把顺序写反了电路照样能编译通过、照样能跑出结果只是结果完全是垃圾。这个坑我踩过一次排查了两个小时才发现是 append 的参数顺序问题。4.3 逆 QFT 与比特序陷阱逆量子傅里叶变换这一块Qiskit 提供了现成的QFT类直接用inverseTrue就行不需要自己搭 H 门加受控相位门的那套线路。但有两个细节必须处理对。第一个是do_swaps。QFT 的教科书实现最后需要一串 SWAP 门把比特顺序翻转回来。Qiskit 的QFT类默认do_swapsTrue会带上这串 SWAP。在模拟器上这没问题但如果你关心线路深度比如后面要上真实硬件可以设成False然后在读测量结果的时候自己把比特串反转。省下的 SWAP 对深度影响不小。第二个是比特序的对应关系。这一点需要说清楚因为它直接影响你怎么把测量出来的字符串翻译成相位。QPE 里第 j 个相位比特施加的是 U^(2^j)j0 对应 U 的一次方jn-1 对应最大的幂次。逆 QFT 之后测量结果的第 j 个比特对应相位二进制小数的第 (n-1-j) 位——也就是说量子比特编号越小对应相位的位权越高不是正好反过来量子比特编号越大对应相位的位权越高。Qiskit 的counts字典里key 是经典寄存器的比特串最左边的字符对应编号最大的经典比特。所以best_key max(counts, keycounts.get) x int(best_key, 2) # 直接按二进制读成整数 phi x / (2 ** n_eval) # 得到 [0,1) 的相位这样读出来的 φ 就是 QPE 估计的相位。我一开始在这里搞混过把比特串反着读结果算出来的能量差了一个 2^n 的因子看起来像是精度不够实际是读法错了。4.4 一份完整可跑的代码把前面所有的东西拼起来这是能直接跑的完整版import numpy as np from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile from qiskit.circuit.library import QFT, PauliEvolutionGate from qiskit.quantum_info import SparsePauliOp from qiskit.synthesis import LieTrotter from qiskit_aer import AerSimulator # ---------- 1. 体系定义 ---------- num_qubits 2 H SparsePauliOp.from_list([ (ZI, -0.4), (IZ, -0.4), (XX, 0.15), (YY, 0.15), (ZZ, 0.1), ]) eigs np.linalg.eigvalsh(H.to_matrix()) E_exact eigs[0] print(f精确基态能量: {E_exact:.6f}) # ---------- 2. 参数推导 ---------- E_ref -1.5 # 保守下界 H_shifted H SparsePauliOp.from_list([(II, -E_ref)]) shifted_eigs np.linalg.eigvalsh(H_shifted.to_matrix()) span shifted_eigs.max() - shifted_eigs.min() t 2 * np.pi / (1.2 * span) # 留 20% 余量 n_eval 10 dE 2 * np.pi / (t * 2 ** n_eval) print(f能量窗口: {2*np.pi/t:.4f}, 分辨率: {dE*1000:.3f} mHa) # ---------- 3. 电路构建 ---------- qr_e QuantumRegister(n_eval, eval) qr_t QuantumRegister(num_qubits, tgt) cr ClassicalRegister(n_eval, c) qc QuantumCircuit(qr_e, qr_t, cr) qc.h(qr_e) for j in range(n_eval): evo PauliEvolutionGate(H_shifted, time-t * (2 ** j), synthesisLieTrotter(reps2)) qc.append(evo.control(1), [qr_e[j]] list(qr_t)) qc.append(QFT(n_eval, inverseTrue, do_swapsTrue).to_gate(), qr_e) qc.measure(qr_e, cr) # ---------- 4. 模拟 ---------- sim AerSimulator() qc_t transpile(qc, sim) result sim.run(qc_t, shots8192).get_counts() # ---------- 5. 结果解析 ---------- best_key max(result, keyresult.get) x int(best_key, 2) phi x / (2 ** n_eval) E_shifted 2 * np.pi * phi / t E_qpe E_shifted E_ref print(f测得相位: {phi:.6f}) print(fQPE 能量: {E_qpe:.6f}) print(f误差: {(E_qpe - E_exact)*1000:.3f} mHa) print(f主峰概率: {result[best_key]/8192:.3f})跑一遍正常情况下误差应该落在 1~3 mHa 这个量级主峰概率在 0.3 到 0.45 之间。如果你的误差是几十个 mHa或者主峰概率低到 0.1 以下那多半是前面说的卷绕问题、比特序问题、或者初态制备有问题。4.5 初态制备决定了你测到的是谁这一段容易被忽略但它直接决定你能不能测到基态。QPE 的前提是输入态是本征态。如果你的初态是多个本征态的叠加测量会以各分量的概率振幅平方为概率随机坍缩到其中一个。所以你能不能拿到基态能量取决于你的初态和基态的重叠度有多大。我上面那段代码里目标寄存器的初态是 |00⟩也就是全部比特置零。这个态和基态的重叠度有多少不知道取决于哈密顿量的具体形式。对某些哈密顿量|00⟩ 和基态的正交性甚至会让你永远测不到基态。标准做法是用 Hartree-Fock 态作为初态。HF 态是占据数意义下能量最低的计算基态通过一组 X 门把对应比特翻转就能制备出来。qiskit-nature的HartreeFock初始化器会自动算出来哪些比特该翻转。from qiskit_nature.second_q.circuit.library import HartreeFock hf HartreeFock(num_spatial_orbitals2, num_particles(1, 1), qubit_mappermapper) qc.compose(hf, qubitslist(qr_t), inplaceTrue)对 H2 这类体系HF 态和真实基态的重叠度通常在 0.9 以上所以主峰概率会相当漂亮。对强关联体系比如拉伸到键长很长的 H2或者某些过渡金属配合物HF 态的重叠度会掉到 0.5 以下这时候你可能需要用多参考态做初态或者接受较低的采样成功率、多跑几轮。提示主峰概率低不一定是电路搭错了。先做个简单的判断——把初态直接设成通过经典对角化得到的基态本征矢如果这样主峰概率能到 0.9 以上而实际初态只有 0.3那就是初态重叠度的问题不是 QPE 的问题。5. 从测量比特串读回能量换算、误差与采样5.1 二进制定点数的读法QPE 输出的相位是 [0,1) 区间内的一个实数用 n 位二进制小数表示。假设测量结果是比特串1011000000n10那么φ 0.1011000000₂ 1/2 1/8 1/16 0.6875换算成能量就是E 2π · 0.6875 / t再减掉 E_ref 的偏移。这里有个必须理解的性质如果真实相位不是 1/2^n 的整数倍QPE 得到的是一个概率分布不是一个确定值。因为相位寄存器只有 n 位精度真实相位落在两个可表示值之间时测量结果会在邻近几个值上分布。经典结论是当真实相位恰好落在两个格点正中间时主峰概率下界是4/π² ≈ 0.405这个数字很重要它解释了为什么你的主峰概率永远上不了 0.5。看到 0.4 左右不要以为是电路有问题这是 QPE 的固有性质。实际处理上我一般不会只取主峰。更好的做法是取概率最高的若干个结果按概率加权求期望或者用中位数 邻域拟合的方式。加权平均能显著降低统计误差尤其在你 shots 有限的时候total sum(counts.values()) E_weights [] for key, cnt in counts.items(): x int(key, 2) phi x / (2 ** n_eval) E 2 * np.pi * phi / t E_ref E_weights.append((E, cnt / total)) E_mean sum(e * w for e, w in E_weights) E_std np.sqrt(sum(w * (e - E_mean) ** 2 for e, w in E_weights)) print(f加权能量: {E_mean:.6f} ± {E_std:.6f})这个做法在相位分布比较集中主峰附近两三个值时非常有效能把有效精度往 2^n 之外再推一点。但如果分布弥散得厉害说明初态重叠度差或者 Trotter 误差大加权平均反而会把误差摊平掩盖问题。5.2 采样次数与主峰概率shots 的选择不是越多越好得看你要什么。如果你只需要主峰对应的比特串argmaxshots 的需求由主峰概率决定。主峰概率 0.4 时要保证主峰在采样结果里排第一用二项分布的粗略估计需要shots ≥ 10 / p_dominantp 0.4 时大约 25 次就够了理论上。但我一般还是跑 4096 到 8192原因是你想要的不只是一个比特串还有整个分布的形状用来做误差估计和异常判断。shots 太少的时候次峰的位置和高度都不可信你没法判断分布是正常的两点弥散还是异常的多峰乱跳。如果你要做加权平均shots 的需求会高一到两个数量级。统计误差大致按1/√shots衰减8192 shots 相比 1024 shots噪声水平能降一个数量级。区分不同场景的推荐值使用场景建议 shots理由快速验证电路是否跑通512只要能看出主峰在哪就行估算单一能量值4096主峰稳定次峰可辨加权平均 误差棒16384分布形状需要足够统计量对比不同 Trotter 步数8192 × 3 组需要能分辨组间差异注意模拟器上的 shots 和真实硬件上的 shots 语义不同。模拟器是精确重复没有硬件噪声所以分布形状反映的是电路本身的数学性质真实硬件上还叠加了门的错误率主峰会被噪声抹平。做方法学研究时用模拟器足够但别把模拟器上的概率分布直接当成硬件上的预期。5.3 Trotter 误差怎么估PauliEvolutionGate内部用的是 Trotter 分解而 Trotter 是近似。误差来自哈密顿量各项之间的不对易性[H_i, H_j] ≠ 0 时e^(-iH_1t)e^(-iH_2t) ≠ e^(-i(H_1H_2)t)一阶 Trotter 的误差量级是 O(t²/reps · Σ|[H_i, H_j]|)二阶是 O(t³/reps²)。这个式子告诉我们两件事演化时间 t 越大误差越大但我们在第 2.3 节已经把 t 固定了所以唯一的调节旋钮就是 reps。实测的做法很直接跑 reps 1、2、4 三组看能量收敛for reps in [1, 2, 4]: qc_test build_qpe_circuit(H_shifted, t, n_eval, num_qubits, reps) E run_and_extract(qc_test) print(freps{reps}: E{E:.6f})如果 reps 从 2 加到 4 能量变化小于 0.5 mHa就可以停了。对小体系、短演化时间通常 reps2 就到收敛区了。对更大的分子或者更长的 t可能需要 4 到 8。还有一个容易被忽略的误差源受控演化门本身的精度。PauliEvolutionGate.control(1)内部会把每个 Pauli 项的指数化电路都加上控制逻辑这个过程中如果某些项分解出的基本门不支持直接受控Qiskit 会做额外的分解可能引入额外的近似。我的经验是如果某个项是单比特 Pauli比如 ZI 这种受控版本很干净如果是多比特 Pauli比如 XX受控后的线路会明显变长这时候 Trotter 的数值误差和线路合成的近似误差会叠加。5.4 比特不够用迭代式 QPE 顶上标准 QPE 有个硬伤为了拿到 2 位精度你得同时保留 n 个相位比特并做逆 QFT电路宽度随精度线性增长还要保证相干性维持在整个过程里。在当前硬件上n 超过 10 就基本没法做了。迭代式 QPEIQPE就是为这个场景设计的。它只用一个相位比特反复测量 n 轮每轮从前一轮的结果里提取相位的一个低位再反馈到下一轮的相位门里。核心思路是第 k 轮测量的是相位的第 k 个二进制位测完之后把已知的高位信息通过一个受控相位门减掉让下一轮能测到更高精度的低位。好处很明确只需要 1 个相位比特 1 个目标寄存器电路宽度从 n m 降到 1 m。代价是电路深度线性增加要跑 n 轮而且每轮都要做一次测量并在线调整下一轮的参数对硬件的实时反馈能力有要求。用 Qiskit 实现的话qiskit-algorithms里的IterativePhaseEstimation可以直接用但如果你要精确控制每一轮的参数尤其是把哈密顿量平移、Trotter 步数这些细节捏在手里自己搭也就几十行。对本地模拟来说标准 QPE 在 n ≤ 12 时完全没问题statevector 模拟的速度完全能接受。只有当你需要 n 16 以上、又想在真实硬件上跑的时候IQPE 才是必要的。6. 踩坑记录与排查速查6.1 常见问题速查表这些是我和身边朋友在跑 QPE 时反复遇到的问题整理成表方便对照排查现象最可能的原因排查方向QPE 能量系统性偏高几十 mHa忘了加核排斥能常数检查nuclear_repulsion_energy有没有加回主峰概率低于 0.15初态与基态重叠度太低换成 Hartree-Fock 初态或检查初态制备门是否加对所有结果都在0...0⟩ 附近受控演化门的时间参数接近 0能量值看起来对但差一个整数倍窗口卷绕能量超出 2π/t 窗口增大 E_ref 的偏移量或减小 t比特串读出的相位和预期差 2^n 倍二进制读法或比特序搞反了先用手算一个已知相位验一下读法增加 reps 后能量反而变差门合成误差超过了 Trotter 误差比较 reps1/2/4 的趋势找拐点模拟跑不动、内存爆掉总比特数超过 20做活性空间约化或者改用 MPS 后端受控演化门 append 后编译报错比特顺序或门的目标比特数不匹配确认.control(1)后是 1 m 个比特6.2 几条只能靠踩坑换来的经验第一先让整条流程在玩具体系上通关再上真实分子。我一开始就是直接拿 H2 的真实积分结果往 QPE 里塞中间任何一环出错都要在十几个环节里排查。后来改成先用一个 2 比特的手写哈密顿量跑通全流程用np.linalg.eigvalsh做绝对参考确认误差在 mHa 量级之后再把哈密顿量换成真实分子的。这样一旦出问题答案就在换什么呢这一步排查范围立刻收敛。第二能量窗口宁可宽一点别抠得太紧。我最初为了追求精度把 t 算得刚刚好让窗口等于能量跨度结果因为活性空间约化后激发态能量估错了那么一点点最高几个本征态直接卷绕出去表现为主峰旁边有个不该出现的次峰。后来改成 1.2 倍余量再没出过这事。多留的那点窗口损失的精度不到 20%但避免的排查时间是几个数量级。第三每次改动参数都存一个基准结果。QPE 涉及的参数太多了——n_eval、t、reps、偏移量、初态、映射方式——改一个就得重跑很容易出现改完 A 发现结果变了但不知道是 A 的影响还是之前 B 就没跑对。我的做法是建一个固定的 reference case每次改动都先跑 reference 确认没退化再跑新配置。这个习惯是从写经典代码时带过来的在这里同样管用。第四别迷信内置的PhaseEstimation类。Qiskit 早期版本里有现成的PhaseEstimation类一行
