4、初始条件与边界条件：初始坐标与速度的生成、周期性边界条件、最小镜像约定、截断半径与邻居列表（Verlet列表、Cell list）

分子动力学（MD）模拟的本质是对牛顿运动方程进行数值积分，而积分过程的起点——初始坐标与初始速度，以及模拟盒子边界的处理方式，直接决定了系统能否快速达到平衡、计算效率的高低以及结果的物理正确性。本章围绕“如何启动一个模拟”和“如何高效且正确地计算相互作用”两大核心问题展开，系统讲解初始条件生成、周期性边界条件（PBC）、最小镜像约定、截断半径以及两种主流邻居列表算法。这些内容构成了MD模拟中“预处理”与“近邻搜索”的基石，是每一位模拟实践者必须精确掌握的工程细节。

上图展示了本章内容的整体工程链路：从初始条件的设定出发，经过周期性边界条件的构建，再到最小镜像约定的应用，最终落实到邻居列表的高效搜索。四个模块环环相扣，任何一环的参数设置不当，都可能导致模拟结果失真或计算效率急剧下降。下面逐一深入剖析。

4.1 初始坐标与速度的生成

4.1.1 核心原理

MD模拟的起点是确定体系中所有原子的位置矢量 r_i(0) 和速度矢量 v_i(0)。初始坐标通常来源于实验晶体结构（如X射线衍射数据）、量子化学优化结构或蒙特卡洛随机堆积。对于液体或无序体系，常用随机放置加能量最小化的策略。初始速度则必须满足麦克斯韦-玻尔兹曼分布，即在温度T下，每个速度分量的概率密度为：

f(v_x) = (m / 2πk_B T)^(1/2) × exp(-m v_x² / 2k_B T)

实际生成时，通常使用Box-Muller变换从均匀随机数生成高斯分布随机数，再乘以目标温度对应的最概然速率。关键约束是：体系总动量必须严格为零，否则系统质心会发生整体漂移，导致温度计算产生偏差。

4.1.2 操作流程与温度标定

1. 坐标预处理：检查原子间是否存在重叠（距离小于0.8Å），若有则进行能量最小化（如最速下降法或共轭梯度法）或软势预平衡。
2. 速度生成：为每个原子生成三个独立的高斯随机数作为速度分量。
3. 动量修正：计算总动量 P_tot = Σ m_iv_i，从每个原子速度中减去 P_tot/M_tot。
4. 温度标定：计算实际温度 T_inst = (2/3Nk_B) × Σ (½m_iv_i²)，若与目标温度偏差超过5%，将所有速度乘以缩放因子 λ = √(T_target/T_inst)。

4.1.3 案例拆解：液态氩体系的初始条件设置

假设我们要模拟一个包含864个氩原子的液态体系，温度85K，密度1.4 g/cm³。具体操作如下：

• 坐标生成：将原子放置在面心立方（FCC）晶格上，晶格常数 a = (4M/ρN_A)^1/3 ≈ 5.36 Å，每个晶胞4个原子，共6×6×6=216个晶胞，恰好864个原子。
• 速度初始化：目标温度85K，氩原子质量 m = 39.948 amu。生成864×3个高斯随机数，修正总动量为零。
• 预平衡：先在NVT系综下运行2000步（时间步长2fs），使用Berendsen恒温器（耦合常数0.1ps），让体系从初始晶格结构“熔化”为液态。
• 结果检查：预平衡后计算径向分布函数g(r)，确认第一峰位置约3.8Å（对应LJ参数σ=3.405Å），配位数约12，符合液态氩特征。

4.2 周期性边界条件与最小镜像约定

4.2.1 核心原理

真实宏观体系包含约10²³个原子，MD模拟只能处理10³~10⁹个原子。为了消除有限尺寸带来的表面效应，引入周期性边界条件（PBC）：模拟盒子在三维空间无限重复，原子i在盒子中的坐标r_i与它在所有镜像盒子中的坐标r_i+nL（n为整数向量）等价。当原子穿过盒子边界时，它从对面边界重新进入，保持原子数守恒。

最小镜像约定是PBC的必然推论：计算原子对i-j的相互作用时，只考虑它们之间距离最近的镜像对。数学上，对每个分量应用：

Δx = x_i - x_j
Δx = Δx - L × round(Δx / L)   // round为四舍五入取整

这样得到的距离向量满足 |Δx| ≤ L/2，三维距离 r = √(Δx²+Δy²+Δz²) ≤ √3·L/2。该约定保证了每个原子对只被计算一次，且计算的是物理上最邻近的相互作用。

4.2.2 判断方法与约束条件

4.2.3 实战场景：纳米液滴的边界处理

模拟一个直径5nm的水纳米液滴（约2000个水分子）在真空中的蒸发行为。此时不能使用3D-PBC，因为周期性复制会导致液滴与自身镜像相互作用，产生虚假的“液滴阵列”。正确做法是：

• 使用一个大盒子（边长≥20nm），仅采用非周期性边界（即自由边界），液滴位于盒子中心。
• 或者使用2D-PBC：在x-y方向周期性，z方向自由，模拟液滴在固体表面的润湿行为。
• 对于需要保持液滴形状的模拟，可施加球形约束势（如软势阱）防止液滴漂移。

4.3 截断半径与邻居列表

4.3.1 截断半径的工程意义

对于Lennard-Jones势，相互作用随距离r增大而快速衰减（~r⁻⁶）。当r超过2.5σ（约8.5Å）时，作用力已小于最大值的1%。因此，设定截断半径r_cut，只计算r < r_cut的原子对，可大幅降低计算量。但直接截断会引入能量不连续，需进行势能平移修正：

U_shifted(r) = U_LJ(r) - U_LJ(r_cut)   (r ≤ r_cut)
U_shifted(r) = 0                        (r > r_cut)

同时，压力计算需加入长程修正项：P_corr = (16/3)πρ²εσ³[(2/3)(σ/r_cut)⁹ - (σ/r_cut)³]。

4.3.2 Verlet列表算法

Verlet列表的核心思想是：为每个原子维护一个“邻居列表”，包含所有距离小于 r_list = r_cut + skin 的原子。由于原子在短时间内移动距离有限，该列表在若干步内保持有效，无需每步重建。

1. 构建阶段：每N_update步（通常10-20步）重建一次列表，遍历所有原子对，计算距离，将r < r_list的原子对加入列表。
2. 使用阶段：在后续N_update步内，只计算列表中的原子对相互作用，跳过距离判断。
3. 参数选择：skin值越大，列表更新频率越低，但列表越长，计算量增加。典型skin = 1.5~2.0Å，N_update = 10~20步。

4.3.3 Cell list（网格列表）算法

Cell list将模拟盒子划分为边长为l_cell ≥ r_cut的网格。每个原子根据其坐标被分配到对应网格单元。计算相互作用时，只需搜索原子所在单元及其相邻的26个单元（三维），而非整个盒子。算法流程：

1. 划分网格：nx = floor(L / l_cell)，l_cell ≥ r_cut
2. 原子分配：cell_id = (ix, iy, iz)，ix = floor(x / l_cell)
3. 邻居搜索：对原子i，遍历其所在cell及26个相邻cell中的原子j
4. 距离判断：仅当r_ij < r_cut时计算相互作用

Cell list的复杂度为O(N)，但需要额外的内存存储网格索引。当体系密度均匀时，每个网格单元原子数约为 ρ×l_cell³，搜索效率极高。

4.3.4 案例拆解：100万原子体系的邻居列表优化

某高分子熔体体系包含1,000,000个原子，盒子边长L=400Å，r_cut=12Å。对比两种方案：

分析：纯Verlet列表构建时需遍历N²/2=5×10¹¹对原子，耗时巨大；Cell list将构建复杂度降为O(N)，但每步重建又增加了距离判断开销。最优方案是Cell list辅助构建Verlet列表：每20步用Cell list快速生成候选对，再构建Verlet列表，后续19步直接使用列表。

4.4 总结复盘与常见误区

本章内容可归纳为“一个中心、两个基本点”：以正确计算相互作用为中心，以初始条件合理性和边界条件有效性为基本点。常见误区包括：

• 误区一：初始速度未修正总动量为零，导致体系质心漂移，温度计算虚高。
• 误区二：盒子边长小于2×r_cut，导致原子与自身镜像相互作用，产生虚假结构。
• 误区三：对库仑力直接截断，忽略长程修正，导致能量误差巨大。
• 误区四：Verlet列表更新间隔过长或skin过小，造成原子对漏算。
• 误区五：Cell list网格尺寸小于r_cut，导致搜索不完整。

正确的工程实践是：在正式模拟前，先进行短时间NVE测试（100-200步），检查能量守恒曲线是否平坦；再进行NVT平衡（10-50ps），监控温度、密度、径向分布函数是否收敛；最后才进行生产模拟。每一步的参数选择都应记录在模拟日志中，确保可重复性。