周期边界下的分子完整性与周期性判断

在分子模拟中，研究体相体系时通常采用三维周期性边界条件，通过有限大小的模拟盒子及其周期镜像，近似描述宏观体系。在这样的坐标表示下，分子的连接可能跨越盒子边界。但仅凭“存在跨边界的连接”，我们还不能判断：它是一个真正的周期性分子，还是一个跨周期边界的有限分子？

这个判断的意义并不限于结构显示。对于正确处理周期相互作用的能量和力计算，只要保持物理连接及相应的镜像关系一致，选择哪一个周期像表示原子或分子，本身不应改变计算结果。分子在存储坐标中被边界分开，也不意味着它在物理上断开。然而，当程序需要计算分子的几何中心、质心，或以这些中心定义拉伸、约束和偏置坐标时，就必须明确哪些原子镜像共同组成同一个分子。直接对被分别映射回盒内的原子坐标取平均，可能得到完全不同的中心；如果这个中心参与施力，错误就会进一步影响模拟。¹

例如，一个蛋白质的部分原子位于盒子的一侧，其余原子通过周期边界与之相连。虽然在当前坐标中它被边界分开，但重新选择各原子的周期镜像后，仍然可以恢复一个完整的有限分子。另一种情况是聚合物通过周期边界首尾连接，沿着化学键可以不断走向更远的周期像，形成无限延伸的链。两者都可能表现出跨边界的键，却不能采用相同的有限分子表示。

因此，对于跨越周期边界、但大小仍然有限的分子，程序需要能够提供一组与连接拓扑一致的完整展开坐标，或者与之等价的镜像位移信息，供质心等分子整体量使用。这里的“完整”是指同一构象中全部连接能够同时得到一致表示；若还要求分子位置随时间连续，则需要进一步跟踪整体镜像的变化。传统分子模拟程序如何提供这种表示，它又在什么情况下不存在？晶体网络的商图方法可以为这个问题给出数学判据。

为什么计算分子中心需要完整坐标

考虑一个边长为 $L$ 的一维周期盒。两个等质量、相互成键的原子被记录在 $0.05L$ 和 $0.95L$。假设这条键连接的是跨边界相距 $0.10L$ 的两个镜像，那么直接平均存储坐标得到的中心 $0.50L$，就落在了盒子中央。将第二个原子表示为 $-0.05L$ 后，这个完整双原子分子的中心应为 $0$，或与之周期等价的 $L$。$0.50L$ 与 $0$ 并不是同一个中心的两个周期像。

对一般有限分子，选定完整展开坐标 $\mathbf R_i$ 后，几何中心和质心分别为

$$\mathbf R_{\mathrm{geom}}=\frac{1}{N}\sum_{i=1}^{N}\mathbf R_i, \qquad \mathbf R_{\mathrm{COM}}=\frac{\sum_{i=1}^{N}m_i\mathbf R_i}{\sum_{i=1}^{N}m_i}.$$

这里 $N$ 是原子数，$m_i$ 是原子质量。若对所有原子施加相同的晶格平移，两个中心也只发生相同的平移；若分别改变各原子的镜像后直接求平均，所得中心却未必与原来的中心周期等价。这说明，计算整体量之前需要先确定分子内部的镜像关系，然后才能讨论整个分子及其中心如何做周期性平移。

这种要求可以通过临时重建坐标实现，也可以通过原子坐标与整数镜像编号共同实现，并不要求动力学积分器始终存储一套展开坐标。关键是调用质心、几何中心等计算的模块，能够使用一致的分子表示。下面从常见的整体映射处理出发，考察这种表示是如何建立的。

传统分子模拟程序中的整体映射与周期性平移

周期性边界条件允许同一个原子用不同晶胞中的镜像来表示。设盒矩阵为 $H=(\mathbf a,\mathbf b,\mathbf c)$，三列分别是晶格矢量，那么 $\mathbf r_i$ 与 $\mathbf r_i+H\mathbf p$（$\mathbf p\in\mathbb Z^3$）表示原子 $i$ 的不同周期像。

对于已经具有完整坐标 $\mathbf R_i$ 的有限分子，整体映射可以通过对全部原子施加相同的晶格平移实现：

$$\mathbf R_i'=\mathbf R_i+H\mathbf k,\qquad \mathbf k\in\mathbb Z^3.$$

统一平移保持分子内部相对坐标不变。$\mathbf k$ 可以根据质心、参考原子或其他分子的位置选择；将分子移回盒内，只是整体镜像选择的一种规则。

但如果存储坐标中的原子分别被映射到了不同晶胞，统一平移整个分子并不能恢复它的内部完整性。此时需要先沿连接关系，为不同原子选择各自的整数位移 $\mathbf p_i$：

$$\mathbf R_i=\mathbf r_i+H\mathbf p_i.$$

沿连接逐步选择原子镜像、恢复完整坐标，通常称为 make whole；按分子进行的 image 则选择整个分子的晶格平移。完整分子不必全部位于一个盒子内，一条长的有限链也可以占据多个盒长。

GROMACS 的 gmx trjconv 中，-pbc whole 重建完整分子，-pbc mol 将分子质心映射到盒内；-pbc nojump 则消除相邻构象之间的跨边界跳跃，需要合适的起始构象。它们分别涉及空间完整性、整体镜像选择和时间连续性。²

AMBER 的 CPPTRAJ 中，image 默认按分子的第一个原子选择镜像，加上 center 后改用质心；autoimage 默认围绕第一个分子、以质心选择镜像，并可调整其他分子的相对放置。内部完整性的修复则可由 fiximagedbonds 沿键逐原子处理。³⁴⁵

模拟运行时，GROMACS 拉伸模块计算普通原子组质心，会为其他原子选择最靠近参考原子的周期像，也可使用前一步质心作为参考位置。¹ 这种参考点约定不等同于沿键重建：对较长的有限分子，它未必保持全部连接；对包含多个分子或片段的原子组，也需要明确整体表示的约定。

这些操作都需要选择原子镜像。但能否重建一个完整的有限分子，还取决于连接本身。

整体映射所依赖的有限分子假设

沿键重建时，可以固定一个根原子，逐步分配相邻原子的镜像。如果连接图是一棵树，路径之间不会冲突；如果有环，绕行一圈后就必须回到原来的镜像位置。若回到的是相邻晶胞中的起点，同一个原子编号便被要求占据两个不同的位置，有限展开因此出现矛盾。

GROMACS 已包含这种检查：mk_mshift() 遍历连接图，mk_grey() 比较不同路径给出的整数镜像位移，并记录 inconsistent shifts。源码将周期性分子列为可能原因，但错误拓扑、异常构象或过长连接也可能导致不一致，因此它并不等于自动完成了周期性分类。⁶ 对无限周期分子的处理仍由用户通过 periodic-molecules = yes 指定，程序不会将其输出为完整的有限分子。⁷

CPPTRAJ 的 fiximagedbonds 则跳过已访问的邻居，没有将闭环残余位移用于周期性分类。因此，得到修复后的坐标并不能证明全部连接都能一致展开；这个结论仅针对该工具。⁵

商图方法与 GROMACS 的位移一致性检查具有相通的数学基础。它进一步将这种条件表述为明确的拓扑判据：在连接标签正确的前提下，判断有限展开是否存在，并确定周期维数。这样，分子坐标的重建与其适用条件就能在同一框架下讨论。

从晶体网络借来一种数学语言

Gao、Wang、Guo 和 Sun 在 2020 年发表的 Determining dimensionalities and multiplicities of crystal nets，系统讨论了如何用商图（quotient graph，QG）识别晶体网络的维数与重数。论文沿用了更早的晶体图论工作，并通过具体结构解释了这种方法的意义。它研究的背景是晶体结构筛选，但其数学语言同样适合描述分子模拟中的周期连接。⁸

在无限周期体系里，每个原子有无穷多个镜像。我们可以把所有平移等价的原子压缩成一个节点，只保留一个晶胞中的原子编号。不过，边不能只记录“原子 $i$ 和原子 $j$ 相连”，还必须记录“连接到 $j$ 的哪一个周期像”。因此，给有向边附上一个整数晶格标签：

$$i\xrightarrow{\mathbf n_{ij}}j, \qquad \mathbf n_{ij}\in\mathbb Z^3.$$

它表示：晶胞 $\mathbf u$ 中的 $i$，连接到晶胞 $\mathbf u+\mathbf n_{ij}$ 中的 $j$。反向走同一条边时，标签取负号。这样得到的是一个有限的、带标签的商图。它可以包含连接同一对节点的多条边，因为这些边通向不同的周期像。

图 1｜石墨烯网络与商图。转载自 Gao 等，2020，原文 Fig. 1，© The Author(s) 2020，依据 CC BY 4.0 使用；图片未作修改，中文说明为本文撰写。原图页面。

图 1a 中的白色和黑色原子分别对应图 1b 的 $n_1$ 和 $n_2$。整个无限石墨烯网络，在这个表示下只需要两个节点与三条边。蓝色边 $e_2$ 位于同一晶胞，标签为 $(0,0)$；沿图中箭头，红色边 $e_1$ 从 $n_1$ 指向 $n_2$，标签为 $(1,0)$，绿色边 $e_3$ 从 $n_2$ 指向 $n_1$，标签为 $(0,1)$。图中的 10、01 表示晶胞编号，不是键长。

沿红边再沿蓝边反向返回，节点编号回到 $n_1$，累计晶格平移却为 $(1,0)$。沿蓝边再沿绿边返回，累计平移为 $(0,1)$。这就是图中的两个基本环 $c_1$ 和 $c_2$。在商图中它们闭合，在真实空间中却分别连接晶胞 00 与 10、00 与 01 中的等价原子。这里的“环”并不是石墨烯中那个熟悉的六元环。

分子的周期性，就是整体展开的可解性

现在可以把最初的问题写成方程。若希望展开后的坐标保留每条边指定的物理连接，就需要为所有原子找到整数向量 $\mathbf p_i$，使每一条边都满足

$$\mathbf p_j-\mathbf p_i=\mathbf n_{ij}.$$

考虑任意闭合路径 $C$。沿它将这些方程相加，左边的每个 $\mathbf p_i$ 都会出现一次正号和一次负号，因此必须有

$$\mathbf w(C)=\sum_{(i,j)\in C}\mathbf n_{ij}=\mathbf 0.$$

反过来，如果所有闭合路径的累计平移都为零，就可以固定一个根原子，沿任意连接路径分配其他原子的 $\mathbf p_i$。两条路径给出的结果之差对应一个闭合路径，其累计平移为零，因而结果与路径无关。这里同时给出了必要性和充分性：

$$\boxed{ \text{存在一致的有限展开} \iff \text{所有闭合路径的晶格平移和为零} }$$

这个结论针对的是连接图中的每个连通分量，并以正确的镜像连接标签为前提。只要存在一个非零的 $\mathbf w(C)$，沿同一路径反复行走，就能到达起点平移 $\mathbf w(C)$、$2\mathbf w(C)$、$3\mathbf w(C)$ 等位置的镜像。对应的无限展开因而包含无穷多个相连原子。

对于同一个连通分量，如果有限展开存在，那么任意两组满足全部边约束的整数位移，只能相差一个对所有原子相同的整数向量。这意味着完整坐标在整体晶格平移之外是确定的，几何中心和质心也随之确定到一个整体周期像。它正好提供了程序所需的性质：原子跨界可以改变存储坐标，却不应任意改变同一个有限分子的内部结构和中心。

由此也能看出，单独一条跨界键为什么不能证明分子具有周期性。一个跨界的水分子可以把两颗氢原子移到氧原子附近；一条没有首尾连接的长链也可以逐段展开。即使普通环状分子跨过了边界，只要沿环一圈的平移相互抵消，仍然是有限分子。真正的区别发生在闭环是否留下非零的晶格平移。

如果还想知道体系沿多少个独立方向延伸，可以收集基本环的平移向量作为矩阵 $S$ 的行，计算

$$d=\operatorname{rank}(S).$$

这里的秩可在有理数或实数上计算。$d=0$ 表示有限分子，$d=1$ 表示周期链，$d=2$ 表示周期层，$d=3$ 表示三维周期网络。图 1 中两个环给出 $(1,0)$ 和 $(0,1)$，所以石墨烯的周期维数为 2。这个维数描述连接向周期像延伸的独立方向数，不是原子坐标是否近似落在一条直线或一个平面上。弯曲的有限链依然是 0D，起伏的周期层也可以是 2D。

为什么复制几个盒子还不够

另一个直观办法是复制出一个超胞，观察连接是否持续向外延伸。但观察范围本身会影响判断：一条连接路径可能需要经过更远的晶胞，才抵达某个等价原子。论文比较了商图方法与依赖有限超胞的既有方法，并用相互穿插的网络说明其中的困难。

图 2｜相互穿插的周期网络。转载自 Gao 等，2020，原文 Fig. 2，© The Author(s) 2020，依据 CC BY 4.0 使用；图片未作修改，中文说明为本文撰写。原图页面。

图 2a 和 2b 分别给出 $\mathrm{Cu_2O}$ 与 $\mathrm{Ag(B(CN)_4)}$ 的结构。图 2c 中的红、蓝两套网络在空间上相互穿插，但彼此没有图中的连接。注意，这里的红蓝代表不同子网；图 2a 中的红蓝则代表氧和铜元素。

对于氧化亚铜，图 2d 的商图给出一组基本环平移：

$$S= \begin{pmatrix} 0&1&1\\ 1&0&1\\ 1&1&0 \end{pmatrix}.$$

它的秩为 3，因此每一套连通子网都是三维周期网络。但任意整数线性组合的三个分量之和都是偶数，无法得到 $(1,0,0)$。也就是说，原点的某个氧原子不能沿这些连接到达它在晶胞 $(1,0,0)$ 中的等价原子；后者属于另一套子网。在这个例子中，$|\det S|=2$ 对应两套平移等价而不相连的子网。

论文指出，若在 $2\times2\times2$ 超胞中直接使用拓扑缩放法的原子数比例，这个体系会给出 $24/6=4$，从而被误判为二维。商图则把两个问题区分开来：秩描述延伸的维数，整数平移生成的子格进一步描述哪些周期像实际可达。计算一般网络的重数时，不能随意挑三条独立向量就取行列式，还需要求出所生成整数子格的基。

对分子模拟而言，这个例子的提醒很直接：空间中靠得近、画面中互相穿过，甚至商图中归属于同一组原子编号，都不等于在无限展开中属于同一个连通分量。

在程序内提供与拓扑一致的分子坐标

把这个判据用于模拟程序，首先要确定“连接”是什么。对于固定拓扑的分子力场，应从模型已有的成键关系出发，并按模型处理替代键的约束。若改用距离阈值或氢键定义边，得到的就是相应接触网络或氢键网络的周期性。普通脂质双层虽然铺满盒子的两个方向，以脂质内部共价键定义连接时，每个脂质仍是有限分子。

其次，需要区分不变的周期连接性质与依赖坐标表示的边标签。如果重新选择原子的坐标代表，令

$$\mathbf r_i'=\mathbf r_i+H\mathbf a_i,$$

那么，为保持同一条物理连接，标签应同步变换为

$$\mathbf n_{ij}'=\mathbf n_{ij}+\mathbf a_i-\mathbf a_j.$$

单条边的标签可以改变，但沿闭合路径相加时，新增项全部抵消，$\mathbf w(C)$ 保持不变。因此，“某条键当前跨不跨盒子”依赖坐标表示，“是否存在非零的闭环平移”才是不随这种表示改变的性质。

可以在固定参考表示下保存标签，也可以从可靠的初始构象推断并跟踪镜像连接。仅有原子编号对的普通 bond list 并不包含全部信息。如果根据坐标推断镜像，则需要保证相应键的镜像选择明确；过长连接、畸变构象和不合适的最近镜像算法，都可能给出错误标签。

对于静态连接图，实现上不必枚举所有环。一次生成树遍历就能分配候选位移 $\mathbf p_i$，随后对每条非树边检查

$$\boldsymbol\Delta_{ij}=\mathbf p_i+\mathbf n_{ij}-\mathbf p_j.$$

全部残余为零，就得到完整展开；存在非零残余，就得到周期绕行的证据。对图的处理复杂度为 $O(N+E)$。若还要区分 1D、2D、3D，则继续求这些残余向量的秩。这里的流程是从上述数学判据导出的实现建议，并非声称 GROMACS 或 AMBER 已按此方式完成周期维数分类。

当有限展开存在时，程序可以把 $\mathbf r_i$ 与 $\mathbf p_i$ 共同作为分子坐标的表示，按需生成 $\mathbf R_i$，提供给质心、几何中心、回转半径或基于这些量的集体变量计算。对整体中心再施加晶格平移，或者对两个中心的距离采用相应的周期边界约定，是之后的步骤，不能代替对分子内部镜像关系的确定。这样，“有限但跨界”的分子仍然能够在程序中作为完整对象参与计算。

若需要计算中心的位移、扩散或随时间连续变化的参考位置，还要跟踪每一步的整体镜像：单帧中完整的分子，在下一帧仍可能整体平移一个盒长。这里需要同时维护空间上的连接一致性和时间上的镜像连续性；仅仅每一步把质心映射回主盒，并不能保证后者。

对真正的周期网络，程序仍然可以选择一个有限片段用于显示，也可以为局部相互作用选择正确镜像；但不能要求每个原子编号只出现一次，同时又在普通欧氏坐标中保留全部周期连接而没有切口。可以定义一个重复单元、局部原子组的中心，或采用专门的周期中心约定，但这些量依赖所选对象和定义，不能直接称为无限网络作为一个有限分子的质心。能否计算某个具体模型的力，还取决于引擎对相应键合项、约束和周期镜像的支持。

因此，区分周期性分子与跨界的有限分子，关系到程序能否为分子整体量提供正确的坐标。对于后者，完整展开使质心等量具有一致的定义；对于前者，则需要保留周期网络的连接性质，并明确整体量所采用的局部或周期约定。闭合路径的晶格平移判据，把这一区分与完整坐标的构造放在了同一个数学框架中。

Footnotes

1. GROMACS，Collective variables: the pull code，尤其是 Definition of the center of mass。文档。查阅日期：2026-09-07。该模块的原子组质心处理说明，镜像约定既用于分析，也用于模拟中的集体变量和施力。 ↩ ↩²

2. GROMACS 2025.1，gmx trjconv，周期性处理选项。文档。 ↩

3. AMBER-hub，CPPTRAJ image。文档。 ↩

4. AMBER-hub，CPPTRAJ autoimage。文档。 ↩

5. Amber-MD/CPPTRAJ，Action_FixImagedBonds.cpp，DoAction()。源码。源码查阅日期：2026-09-07；此处讨论的是轨迹处理工具，不据此推断全部 AMBER 动力学引擎的行为。 ↩ ↩²

6. GROMACS，src/gromacs/pbcutil/mshift.cpp 中的 mk_mshift()、mk_grey() 和 shift_x()。源码。源码查阅日期：2026-09-07；此处用于说明分子重建的路径一致性检查，不代表所有执行路径都使用同一套代码。 ↩

7. GROMACS，Molecular dynamics parameters，periodic-molecules。文档。查阅日期：2026-09-07。 ↩

8. Gao H, Wang J, Guo Z, Sun J. Determining dimensionalities and multiplicities of crystal nets. npj Computational Materials. 2020;6:143. doi:10.1038/s41524-020-00409-0。商图与维数见原文公式 (2)–(6)，氧化亚铜和重数讨论见公式 (7)–(8) 及 Fig. 2。 ↩