周期性边界条件

定义

周期性边界条件是给模拟体系规定的一种环绕式约束：以一块有限体积的模拟盒子作为基本单元，让人为复制出的众多同形副本填满整个空间，这些副本被称为镜像。^[1] 参与计算的粒子实际上代表一个无限周期体系中的一小块体积，盒子内部的粒子既与盒内其他粒子作用，也与各镜像中的粒子作用。^[2]

在固体物理和电子结构计算中，同一思想通常称为玻恩－冯卡门边界条件（Born–von Kármán boundary condition），也常被叫作循环边界条件，其内容是要求波函数在布拉维格点上具有周期性，用于模拟理想晶体。^[3][4] 对晶体点阵振动的处理则要求原子位移以晶体的长度 L = Na 为周期，相当于把晶体首尾相接成一个环。^[5]

原理

在这种边界条件下，主盒子里的粒子移动时，每个镜像盒中的对应粒子以完全相同的方式移动；一旦粒子从某一面穿出盒子，它的一个副本就会从对面的相应位置进入。^[2] 模拟过程中只需保存主盒子内粒子的坐标，其余副本的位置由平移关系推得，因此体系不存在壁面，粒子数密度在模拟中守恒。^[6] 从拓扑上看，二维周期体系相当于把平面卷成一个环面。^[1]

由于副本的存在，两个粒子之间的路径不再唯一，距离也就有多种取法。为此引入最小镜像约定：中心粒子只与最靠近它的那一个镜像发生相互作用，从而把无限多项求和简化为盒内 N-1 个最近邻居的求和。^[7] 对边长为 L 的立方盒子，处理办法是把每个方向的坐标差减去盒长的整数倍，使其落在 $-L/2$ 与 $L/2$ 之间。^[2]

要让最近镜像的算法不出错，截断半径与盒子尺寸必须匹配。如果截断半径超过盒子最小边长的一半，同一个粒子的两个不同镜像可能同时落入截断球内而被重复计算。^[8] 实践中通常要求盒子的各边尺寸不小于截断半径的两倍，而截断半径本身多取 1 纳米左右。^[2]

库仑力这类长程作用不能只保留最近镜像，因为盒内粒子与远处无数镜像之间的微小贡献累加起来不可忽略。^[9] 1921 年，埃瓦尔德提出一种改写办法，把 $1/r$ 势的求和拆成实空间与倒空间两部分，两者都迅速收敛，可以准确计入全部镜像的贡献。^[10] 这一方法被称为埃瓦尔德求和，原先用于晶体点阵势的计算，后来成为周期性体系静电处理的通用手段，并衍生出粒子网格埃瓦尔德（PME）等快速算法。^[11]

在电子结构计算中，周期性边界条件要求晶体势能满足晶格平移对称性，单电子波函数由此满足布洛赫定理，可以写成平面波与周期函数之积 $\psi_{n\mathbf{k}}(\mathbf{r}) = e^{i\mathbf{k}\cdot\mathbf{r}}u_{n\mathbf{k}}(\mathbf{r})$ 的形式。^[3] 波矢相差一个倒格矢的两个状态彼此等价，因此只需在第一布里渊区内取有限的一组 k 点做加权求和，就能近似整个布里渊区的积分。^[12]

发展历程

这种边界条件最早出现在固体物理中。1912 年，马克斯·玻恩与西奥多·冯·卡门在研究晶体点阵振动时首次采用它，把有限长的晶体链处理为无限长链的一个重复单元。^[5] 他们在 1912 年与 1913 年发表的一系列论文中给出了最早的、以晶体假设为基础的固体比热理论之一，这一条件就包含在其中。^[3]

引入该条件后，格波的波矢不再连续，只能取一系列分立值，相邻两个取值之间相差 $2\pi/(Na)$，其中 N 为晶体链的原胞数，a 为晶格常数；第一布里渊区内波矢的取值总数恰好等于原胞数。^[5]

玻恩与冯卡门的这项工作起初未受重视，原因是同年发表的德拜热容理论采用弹性波近似，已经可以较好解释当时的实验；此后更精确的测量显示德拜模型的不足，1935 年布莱克曼重新采用玻恩－冯卡门近似讨论晶格振动，逐步发展成今天的晶格动力学理论。^[5]

计算机模拟兴起后，这一思想被移植到原子尺度的计算中。分子动力学方法在 20 世纪 50 年代末被用于研究硬球体系，1953 年 Metropolis 等人在蒙特卡罗计算中引入最小镜像约定以减少运算量，此后周期性边界条件连同该约定成为分子模拟的常规设置。^[13][2]

应用

分子动力学与蒙特卡罗方法用这一边界条件计算体相气体、液体、晶体和混合物的性质，也常用于研究处在显式溶剂环境中的生物大分子。^[1] 在研究纳米材料的力学行为时，例如对石墨烯做拉伸模拟，常在平面内方向施加周期性边界条件以消除自由边缘的影响，让只有几百个原子的模型表现出无限大材料的行为。^[2]

在电子结构计算中，周期性边界条件是平面波基组方法的前提，晶体、表面与低维材料的计算都建立在这一条件之上。表面体系一般用板层模型处理，即在两个方向上保持周期、在第三个方向留出足够厚的真空层。^[13] 计算时需要确定 k 点网格的疏密与平面波截断能，二者共同决定结果的精度和计算量。^[12]

只在部分方向需要周期的体系可以使用半周期边界条件。例如在 NAMD 软件中模拟金属有机框架 UiO-66 时，x、y 方向取周期、z 方向保持自由，全周期体系用粒子网格埃瓦尔德方法处理静电，半周期体系则改用多层求和法。^[14]

模拟盒子除立方体外，还可取长方体、截角八面体、六方棱柱或三斜盒子等形式，但非立方盒子的最小镜像处理方式与立方盒子不同，需要专门实现。^[15] 研究剪切流动时还用到 Lees–Edwards 边界条件，此时上下镜像盒子以恒定速度相对平移；用于非平衡现象的周期壳边界条件则通过在外围设置镜像粒子来减弱有限模拟区域的畸变。^[16]

局限

人为周期性的引入会留下痕迹。以无序液体为例，宏观上被强加了一个盒子尺度的规则结构，径向分布函数中可能出现并不真实的波纹，在距离接近盒子半边长的位置尤为明显；改变盒子尺寸重复模拟，可以把真实的局域结构与人造特征区分开来。^[9]

有限盒子带来的有限尺寸效应也无法回避。当盒子边长为 L 时，波长大于 L 的密度涨落不会出现，特征尺度超过盒子尺寸的过程（例如某些相变）难以被正确描述。^[6] 若体系受到外场作用，而该外场的势函数不具有与盒子相同的周期性，则这一边界条件不再适用。^[2]

长程静电作用需要额外处理。对净电荷不为零的周期体系，总能量会随所取镜像数目的增加而增大，在无限镜像极限下不具有意义。^[11] 在路径积分分子动力学中，施加这一条件还要考虑聚合物环跨越盒子边界的缠绕构型，忽略缠绕或处理不当会引入偏差。^[17]

体系尺寸对结果的影响在软物质中尤其显著。对膜状体系而言，模拟面积的大小会改变所得性质，例如脂质双层的侧向扩散系数随体系尺寸不同而变化，因此选取合适的模拟尺度本身就是一个需要权衡的问题。^[18] 为减弱立方晶格形状带来的伪影，研究者也提出了球状边界条件等替代方案，用不带晶格结构的边界包裹模拟区域。^[19]

参见

• 最小镜像约定 —— 周期性体系中确定粒子间距离与相互作用对象的常用约定。

• 埃瓦尔德求和 —— 精确处理周期性体系中长程库仑相互作用的经典方法。

• 布洛赫定理 —— 周期性势场中电子波函数所满足的基本定理，是周期性电子结构计算的基础。

• 分子动力学 —— 最常使用周期性边界条件的模拟方法之一。

• 蒙特卡罗方法 —— 另一种同样常与周期性边界条件配合使用的模拟方法。

• 有限尺寸效应 —— 由模拟盒子尺寸有限而引起的系统偏差。