水分子动力学模拟

从分子到电路

问题与输入

先从一个水分子看起：一个氧原子、两个氢原子，各自占据一个位置。如果把三个原子的位置连续记录下来，就得到一条分子轨迹。不过，要从当前的位置算出下一步，还需要知道它们受到怎样的力。

这里的计算从能量入手：用训练好的量子-经典代理势能模型预测能量，再由能量的变化求力，最后更新原子位置。理解这一点，也就能区分本教程与严格的“从头算分子动力学”（AIMD）：后者通常在每一步用电子结构方法计算能量和力，本教程则使用代理模型，并不在每一步重新求解电子结构。

分子构型: 用每个原子的 x、y、z 坐标，就能描述此刻的构型。要看形状怎样变化，可以先看两条 O–H 键长和 H–O–H 夹角：它们反映伸缩与弯曲，也不会因分子的整体平移或转动而改变。
势能与力: 每组原子位置都对应一个能量值 E，这种对应关系构成势能面。沿位置的变化观察能量，便能由负梯度得到力（F = −∇E）；积分器再结合力与原子的运动状态，计算下一时间步的位置。
量子-经典模型: 为了得到这个能量，先把分子几何变成电路输入，再把量子电路的读出整理成数值特征，交给训练好的经典模型预测势能。量子电路在这里承担特征计算，原子位置则在经典计算阶段更新。

用内部几何描述水分子的形状

三个原子的坐标已经给出，怎样让模型关注水分子本身的形状？先从坐标中取出两条 O–H 键长和 H–O–H 夹角，再把它们组合成三个几何特征：

q = (r_{1} + r_{2}, {(r_{1} - r_{2})}^{2}, \cos θ)

这里，r₁、r₂ 是以 Å 为单位的两条 O–H 键长，θ 是 H–O–H 夹角。顺着公式看，三项依次描述总键长、键长差异和弯曲程度。即使整体平移、转动分子，或交换两个氢原子，这三个特征也不会改变。

把几何特征转换为三个旋转角

形状有了数值描述，还需要把它转换成量子门能接收的角度。沿用训练时的均值和尺度，先将几何特征标准化，再做缩放与平移，就得到电路最前面三个 Ry 门的输入角：

φ_{j} = b_{j} + a_{j} \frac{q_{j} - μ_{j}}{s_{j}} (j = 1, 2, 3)

公式中的 μ、s 是训练后固定的几何均值与尺度。默认偏移 b 为 π/2，三项缩放 a 依次为 π/4、π/8、π/4。代入当前构型的特征，就会得到以弧度表示的 φ，也就是下方代码中的 enc。

这里有两种角度需要分清：θ 描述水分子的键角，φ 则是送入电路的旋转角。后续量子门另有 11 个模型参数，对应代码的 theta；读代码时，可以将它们与随构型变化的编码角分别理解。

量子电路提取特征，经典模型预测势能

编码角进入三比特电路之后，离势能还差一步。电路先对选定的 Pauli 算符求期望值，得到量子特征 z；这些特征经过标准化，再交给经典模型，才得到势能预测：

E (R) = μ_{E} + s_{E} f_{w} (\frac{z (R) - μ_{z}}{s_{z}})

R 表示三个原子的全部坐标。对于这个构型，默认模型读出 7 个 Z 基和 7 个 X 基特征（共 14 个），交给经典神经网络 f_w。这个网络含两个 32 单元隐藏层，使用 SiLU 激活。

公式中，μ_z、s_z 先按分量标准化特征，μ_E、s_E 再把网络输出还原为以 eV 为单位的能量。这些参数在模型训练完成后保持固定。还需要记住这里的预测对象：默认训练目标是数据集中的相对势能，量子读出本身不是电子能量。

比较相邻构型的能量，求出原子受力

现在能计算一个构型的势能了，怎样从中得到力？关键是看位置发生微小变化时，能量怎样改变：力就是势能对坐标的负梯度。当前水分子动力学模拟程序分别对每个原子的 x、y、z 坐标施加正、负扰动，比较两侧的能量，用中心有限差分近似求出受力：

F_{i α}^{FD} = - \frac{E (R + h e_{i α}) - E (R - h e_{i α})}{2 h}

读公式时，i 标记原子，α 表示坐标方向，e_iα 只在该坐标上取 1。取默认扰动步长 h = 0.001 Å，将两侧能量之差除以两倍步长，就近似得到了这个方向上的能量变化率；再取负号，便得到以 eV/Å 为单位的力。

力算出后，代码还会移除平移与转动方向的数值残差，再把处理后的力交给 VelocityVerlet 积分器。由它更新位置与速度，便从当前构型走到了下一步。

公式对应的代码与默认配置

把公式与实现对照起来看，可以从下面的代码入口逐项查阅。这里讲解的是单水分子应用的 F2/A2 默认模型；要运行这个模型，还需要从配套 checkpoint 加载模型权重与归一化状态。

几何、编码与量子读出：water_symmetric_invariants、water_symmetric_angle_features、feature_tensor_from_angles。
经典能量模型：_predict_tensor、_fit_normalizers。
坐标差分与刚体残差处理：calculate_geometry_energy_and_force。
默认配置：quantum.encoding、classical、force。

分子构型3 个编码角3 比特电路经典能量模型电路另有 11 个模型参数

带着这条计算路径往下读，代码 1、2 会先搭好可复用的三比特电路。到代码 3，再用当前仓库的默认模型，在 CPU 上算出一个构型的能量与力，把上面的公式落实为数值。

最后的代码 4 会转向另一份归档的轨迹 CSV。它没有附上配套的模型权重、编码配置或执行后端记录，所以目前无法确认它是否由当前模型生成，也无法用前三段代码重现后面的能量曲线。阅读时，需要把这份轨迹与前面的模型计算分开理解。

PivotQ / Python

构造量子电路

几何已经可以写成三个编码角，接下来要把它们放进怎样的电路？先准备 Python 3.12，在 PivotQ 仓库根目录运行 python -m pip install ./packages/framework，也可以使用仓库已配置的统一环境。PivotQ 的电路接口直接复用 Qiskit 对象。

按代码 1、代码 2 的顺序运行，就会先定义辅助函数，再构造三比特电路并打印结构摘要。下方的“代码 2 输出”来自这两段源码的实际运行；浏览器只展示源码和输出，不执行 Python。

此时搭好的是电路结构。它保留了 enc 和 theta 共 14 个未赋值参数，其中 enc 随分子构型改变，theta 是模型训练后使用的参数。摘要中的基础门总数 70、电路深度 39，描述的也都是电路结构，不能表示某一帧构型的能量。

看电路图时，可以从左边沿着 q0、q1、q2 三条水平线往右读。单比特旋转先改变各量子比特的状态，CZ 门再对相邻的两个比特施加受控相位变换；随后，一组 Pauli 旋转在指定比特上沿指定轴继续变换量子态。

两张图展示了同一电路的不同层次：“模型门结构”保留 Rx 与 Pauli 旋转，便于理解；“基础门分解”把它们展开为 Ry、Rz、CZ 门。要让这套结构给出能量，还需要代入参数、读出量子态，再将读出特征交给经典模型。

代码 1导入与辅助函数展开代码

import math
from uuid import NAMESPACE_URL, uuid5

from pivotq import QuantumCircuit, Parameter

ADAPT_OPERATORS = ("IYZ", "YII", "YZI", "IIX", "YII")

def _parameter(name: str) -> Parameter:
    """Create a stable parameter identity so repeated exports are reproducible."""
    return Parameter(name, uuid=uuid5(NAMESPACE_URL, f"qhai-h2o-f2a2/{name}"))


def _native_h(circuit: QuantumCircuit, qubit: int) -> None:
    circuit.rz(math.pi, qubit)
    circuit.ry(math.pi / 2.0, qubit)


def _native_h_inverse(circuit: QuantumCircuit, qubit: int) -> None:
    circuit.ry(-math.pi / 2.0, qubit)
    circuit.rz(-math.pi, qubit)


def _native_rx(circuit: QuantumCircuit, angle, qubit: int) -> None:
    circuit.rz(math.pi / 2.0, qubit)
    circuit.ry(angle, qubit)
    circuit.rz(-math.pi / 2.0, qubit)


def _native_rx_inverse(circuit: QuantumCircuit, angle, qubit: int) -> None:
    circuit.rz(math.pi / 2.0, qubit)
    circuit.ry(-angle, qubit)
    circuit.rz(-math.pi / 2.0, qubit)


def _native_cnot(circuit: QuantumCircuit, control: int, target: int) -> None:
    _native_h(circuit, target)
    circuit.cz(control, target)
    _native_h(circuit, target)


def _native_cnot_inverse(
    circuit: QuantumCircuit, control: int, target: int
) -> None:
    _native_h_inverse(circuit, target)
    circuit.cz(control, target)
    _native_h_inverse(circuit, target)


def _pauli_rotation(circuit: QuantumCircuit, word: str, angle) -> None:
    """Compile exp(-i angle P/2) into the frozen Ry/Rz/CZ gate set."""
    support = [qubit for qubit, symbol in enumerate(word) if symbol != "I"]
    if not support:
        raise ValueError("Identity is not a parameterized operator")
    if len(support) == 2 and support == [0, 2]:
        raise ValueError("Direct q0-q2 Pauli rotations are not allowed")

    for qubit in support:
        symbol = word[qubit]
        if symbol == "X":
            _native_h(circuit, qubit)
        elif symbol == "Y":
            _native_rx(circuit, math.pi / 2.0, qubit)
        elif symbol != "Z":
            raise ValueError(f"Unsupported Pauli symbol: {symbol}")

    for first, second in zip(support[:-1], support[1:]):
        _native_cnot(circuit, first, second)
    circuit.rz(angle, support[-1])
    for first, second in reversed(list(zip(support[:-1], support[1:]))):
        _native_cnot_inverse(circuit, first, second)

    for qubit in reversed(support):
        symbol = word[qubit]
        if symbol == "X":
            _native_h_inverse(circuit, qubit)
        elif symbol == "Y":
            _native_rx_inverse(circuit, math.pi / 2.0, qubit)

代码 2

构建三比特电路

Python

def build_circuit() -> QuantumCircuit:
    """Return the project's frozen, unbound three-qubit logical circuit."""
    circuit = QuantumCircuit(3, name="h2o_f2a2")
    encoding = tuple(_parameter(f"enc[{index}]") for index in range(3))
    theta = tuple(_parameter(f"theta[{index}]") for index in range(11))

    # Geometry encoding: angles from angles.csv map one-to-one to q0, q1, q2.
    for qubit, angle in enumerate(encoding):
        circuit.ry(angle, qubit)

    # Native seed: three Ry rotations, linear CZ connectivity, three native Rx.
    for qubit in range(3):
        circuit.ry(theta[qubit], qubit)
    circuit.cz(0, 1)
    circuit.cz(1, 2)
    for qubit in range(3):
        _native_rx(circuit, theta[3 + qubit], qubit)

    # Frozen ADAPT sequence: IYZ -> YII -> YZI -> IIX -> YII.
    for index, word in enumerate(ADAPT_OPERATORS):
        _pauli_rotation(circuit, word, theta[6 + index])

    circuit.metadata = {
        "model": "H2O_F2A2_R1_DROP_XXX_MLP_ONLY",
        "logical_qubit_order": (
            "q0,q1,q2 correspond to Pauli word characters left-to-right"
        ),
        "adapt_operators": list(ADAPT_OPERATORS),
    }
    return circuit

circuit = build_circuit()
print(f"量子比特: {circuit.num_qubits}")
print(f"未赋值参数: {len(circuit.parameters)}")
print(f"基础门总数: {circuit.size()}")
print("门计数: " + ", ".join(
    f"{name.upper()}={count}" for name, count in sorted(circuit.count_ops().items())
))
print(f"电路深度: {circuit.depth()}")

代码 2 输出 / PivotQ 电路

量子比特: 3
未赋值参数: 14
基础门总数: 70
门计数: CZ=6, RY=27, RZ=37
电路深度: 39

读到深度 39 时，可以把它理解为量子门按依赖关系排成的层数。它描述电路，而分子动力学的时间步描述分子演化，因此这 39 层不能换算成后文的 fs。

图 1

模型门结构

3 个几何编码 Ry 门 → 3 个参数化 Ry 门 → 2 个 CZ 门 → 3 个 Rx 门 → 5 个 Pauli 旋转。

图 2

基础门分解

Ry(角度): 这个门旋转一个量子比特。回看电路的最前面，三个 Ry 分别接收从分子几何得到的 enc[0]、enc[1]、enc[2]。
CZ: 图中连接相邻两个量子比特的是 CZ 门。当两者都处在 |1⟩ 分量时，它会改变该分量的相位。
Pauli 旋转: 读这些字母，就能知道每个比特的作用轴。例如 IYZ 中，I 表示 q0 不参与，Y、Z 分别作用于 q1、q2。
基础门分解: 把逻辑 Rx 和 Pauli 旋转展开为 Ry、Rz、CZ 门，就得到图 2。门的表达变了，仍然是图 1 的同一模型。

CPU + QPU

计算步骤

电路准备好以后，它怎样参与一次位置更新？先代入编码角和模型参数，量子阶段就可以根据 Z 基、X 基的读出结果，各提取 7 个统计特征。这里得到的还只是特征，并不是原子坐标或受力。

这些特征随后进入原应用的经典能量模型。模型有两层、每层 32 个单元，负责预测势能。有了能量计算方法，就能对坐标做中心有限差分，求力并移除刚体数值残差；VelocityVerlet 积分器再结合所得的力，更新原子位置和速度。新构型由此成为下一时间步的输入，同一条计算链便可以继续。

沿着下面的流程图看，QPU 是量子特征阶段选择的目标设备，CPU 则承担经典能量推理、求力与位置更新。

一个时间步内的计算如何衔接

01
原子位置 → 编码角
从当前氧、氢原子的坐标提取分子几何，再沿用模型训练时的编码规则，得到三个 Ry 门角度。后面回放的归档轨迹没有附上每帧的编码角。
02
电路 → 量子特征
把这些编码角连同模型参数代入电路，分别在 Z、X 基读出，就能各取得 7 个供经典模型使用的特征。
03
特征 → 能量与力
经典模型先预测当前构型的势能。再对原子坐标做小幅扰动、比较能量变化，便能通过差分近似求出各原子的力。
04
力 → 下一步位置
力算出之后，VelocityVerlet 积分器结合原子质量及当前速度更新位置。对这个新构型重复计算，就进入下一步。

回头看代码 2，它打印的比特数、参数数和门数，描述的是第 2 步将用到的电路。那段代码还没有代入一帧的参数，也没有进行读出、能量预测或轨迹积分。

读到这里，还可以停下来分清“测量结果”和“读出特征”。一次测量只给出一组比特结果；这里需要的特征，则是根据电路在 Z 基、X 基下的读出结果计算的统计量。在采样设备上，要通过多次测量估计这些统计量；理想模拟器也可以直接计算期望值。

因此，一次测量的结果还不能直接作为分子的势能。要等经典能量模型接收这些特征，才会得到当前构型的能量预测。

执行流程

01
初始化与轨迹控制
目标 CPU
02
量子特征计算
目标 QPU
03
经典能量预测
目标 CPU
04
求力与轨迹积分
目标 CPU
05
轨迹分析
目标 CPU

代码 3

在 CPU 上计算能量与受力

现在把这条计算链落实到一个构型上。使用仓库已配置的 Python 环境，在 applications/h2o-hybrid-aimd/ 目录运行下面这段代码：它会读取默认配置与已训练模型，算出初始水分子构型的势能和三个原子的受力。下方保留了实际运行输出，浏览器不执行 Python。

Python · CPU

import torch
from single_h20_aimd import load_hybrid_potential
from single_h20_aimd.configuration import load_config, project_path
from single_h20_aimd.data import water_internal_to_cartesian

config = load_config("configs/h2o_aimd.yaml")
checkpoint = project_path(config, config["checkpoint"]["path"])
potential = load_hybrid_potential(config, checkpoint)
geometry = water_internal_to_cartesian(*config["aimd"]["initial_internal_coordinates"])

with torch.no_grad():  # MLP 推理 → 中心差分求力 → 刚体残差投影
    result = potential.predict_geometry_energy_and_force(geometry)

print(f"势能: {result.energies_eV[0]:.6f} eV")
print("受力 (eV/Å)，依次为 O、H、H 的 x、y、z 分量：")
print(result.forces_eV_per_A[0].round(6))

运行输出 / CPU · 单个构型

势能: 0.327675 eV
受力 (eV/Å)，依次为 O、H、H 的 x、y、z 分量：
[[-4.251411  0.       -3.290627]
 [ 0.711873  0.        3.471972]
 [ 3.539538 -0.       -0.181345]]

这次计算中，CPU 负责 MLP 能量推理、中心差分求力和刚体残差投影。为了求出力，程序以 0.001 Å 为差分步长，对 9 个坐标分别施加正、负扰动。这样，连同中心构型，一共需要计算 19 个构型的势能。

到这里，已经得到了一次能量与力的结果。原子的位置尚未在这段代码中推进，后续的位置更新由积分器完成。

已保存的数据

轨迹与能量曲线

如果把每一步的位置连起来，会看到怎样的运动？为了练习读轨迹，这里换用另一份归档的 CSV，它记录了每一帧的原子坐标、能量和时间。这份数据并非代码 3 的运行结果；下方的分子视图和曲线只是读取归档，并同步显示同一帧。

先看时间轴：它记录了 0–100 fs 的演化过程。fs 是飞秒，1 fs 等于 10⁻¹⁵ 秒。数据每隔 0.1 fs 保存一帧，因此 1000 个时间步连同初始构型，共有 1001 帧。

再把目光移到分子形状。两条 O–H 键长的增减对应伸缩，H–O–H 夹角的变化对应弯曲。键长用 Å（埃，1 Å = 0.1 nm）表示，夹角用度表示。沿着时间轴看这三个量，就能跟上水分子内部几何的变化。

与结构相对应，能量图用 eV（电子伏特）记录势能、动能和总能量：势能对应当前原子排布，动能对应原子的运动，总能量是两者之和；温度曲线也反映运动强弱。把同一时刻的分子构型与曲线对齐，就可以观察结构变化怎样与能量、温度变化相伴出现。

读一帧：先对齐结构与曲线

先停在第 0 帧：两条 O–H 键都是 0.9572 Å，H–O–H 夹角为 104.52°，温度为 300 K，总能量为 0.5912 eV。这些数值共同描述了轨迹的起点。

把曲线光标拖到 100 fs，两条键长变为 1.2024 Å 和 1.1422 Å，夹角为 95.27°。与起点相比，两条键分别长了 0.2452 Å、0.1850 Å，夹角小了 9.25°。此时，分子图显示当前形状，曲线光标也指向同一帧。

这次对比只告诉了首末两帧的差异。要知道两条键在中间是否一直伸长或收缩，还需要沿时间轴查看完整过程。

怎样看分子形变？

观察时，可以先分开看伸缩与弯曲：两条 O–H 键可能一起伸缩，也可能一长一短；H–O–H 夹角的改变则对应弯曲。实际轨迹会把这些变化叠在一起，所以只看首末两帧，还不能判断中途的振动模式或频率。

总能量应当保持不变吗？

要回答这个问题，需要先知道模拟条件。在理想的封闭、无恒温器模拟中，总能量应大体保持不变；恒温控制与数值积分都可能改变曲线。这份归档缺少原始系综和积分器设置，因此目前还不能凭这条曲线判断能量守恒或计算精度。

代码 4

根据坐标复算首末两帧的键长和键角

图中的键长和键角，也可以直接从原子坐标复算。先在下方的“数据来源与运行条件”下载两个 CSV，放到同一个 imported/ 目录；再把这段代码保存为 analyze_imported.py，在 Python 3.9+ 中运行 python3 analyze_imported.py imported。它会从坐标算出键长和键角，并读取相同时间步的总能量，便于与图中数据对照。

Python · 标准库

"""Inspect the first and last frames of the imported H2O trajectory."""

import csv
import math
import sys
from pathlib import Path


data_dir = Path(sys.argv[1]) if len(sys.argv) > 1 else Path("imported")


def rows(filename):
    with (data_dir / filename).open(newline="", encoding="utf-8") as handle:
        return list(csv.DictReader(handle))


def atom(position, name):
    return [float(position[f"{name}_{axis}_A"]) for axis in "xyz"]


logs = rows("md_log.csv")
positions = rows("positions.csv")
if not logs or len(logs) != len(positions):
    raise ValueError("Energy log and position file must have the same nonzero frame count")

for index in (0, -1):
    log, position = logs[index], positions[index]
    if log["step"] != position["step"]:
        raise ValueError("Energy log and position steps do not match")
    oxygen, hydrogen_1, hydrogen_2 = (atom(position, name) for name in ("O", "H1", "H2"))
    bond_1 = [h - o for h, o in zip(hydrogen_1, oxygen)]
    bond_2 = [h - o for h, o in zip(hydrogen_2, oxygen)]
    length_1, length_2 = math.dist(oxygen, hydrogen_1), math.dist(oxygen, hydrogen_2)
    cosine = sum(a * b for a, b in zip(bond_1, bond_2)) / (length_1 * length_2)
    angle = math.degrees(math.acos(max(-1.0, min(1.0, cosine))))
    print(f"step={log['step']} time_fs={float(log['time_fs']):.1f}")
    print(f"  oh1_A={length_1:.4f} oh2_A={length_2:.4f} hoh_deg={angle:.2f}")
    print(f"  total_eV={float(log['total_energy_eV']):.4f}")

change = float(logs[-1]["total_energy_eV"]) - float(logs[0]["total_energy_eV"])
print(f"delta_total_eV={change:.4f}")

运行输出 / 导入轨迹

step=0 time_fs=0.0
  oh1_A=0.9572 oh2_A=0.9572 hoh_deg=104.52
  total_eV=0.5912
step=1000 time_fs=100.0
  oh1_A=1.2024 oh2_A=1.1422 hoh_deg=95.27
  total_eV=0.4855
delta_total_eV=-0.1056

从输出可以读到，首末帧的总能量相差 −0.1056 eV。这个差值是一个观察结果；要进一步判断能量是否守恒，还需要查看整条时间序列和模拟条件。现有归档也没有原始运行后端记录，因此无法据此判断变化来自物理过程、数值方法还是硬件。

查看水分子的运动

数据来源与运行条件

回放所用的是已导入的教程轨迹：1000 步、100 fs，共 1001 帧，初始温度为 300 K。导入文件没有附上原始运行编号与计算后端，因而这些轨迹可以用于学习如何读数据，但不能用于比较硬件性能。

追溯数据时，还需要留意：这组轨迹、工作台教程中的 CPU 操作教程以及另行保留的历史运行记录，分别存档。具体来源可以继续查看数据来源说明。

下载能量与时间数据 CSV 下载原子坐标 CSV

.wb[data-aimd-workbench] > .wb-load-status, .wb[data-aimd-workbench] > .wb-analysis-grid { display: none; }曲线和三维轨迹需要 JavaScript。下载能量与时间数据 CSV 下载原子坐标 CSV

正在读取轨迹数据…

水分子三维轨迹

OH₁H₂

显示运动轨迹

O–H₁—Å

O–H₂—Å

H–O–H—°

动力学曲线

时间范围