Long-time Brownian dynamics of hard and sticky colloids by geometric projection

技术摘要：通过几何投影实现硬球与粘性胶体的高效布朗动力学模拟

问题陈述
使用布朗动力学（BD）进行胶体悬浮液的粒子级模拟经常受到数值模型中极快微观时间尺度的限制。虽然标准布朗动力学消除了惯性弛豫，但通常仍保留显式的近接触势（陡峭的排斥壁和狭窄的吸引阱）来模拟硬球排斥和粘性相互作用。这些势能引入了快速的位置弛豫时间，迫使时间积分器使用极小的时间步长。这种限制对于聚合胶体（其中不可逆的粘附接触会持续多个步骤）以及高密度系统（其中存在多体重叠事件）而言尤为严重。因此，诸如聚合动力学、凝胶化和电泳沉积等感兴趣的集体过程，由于难以承受解析这些快速近接触动力学的巨大计算成本，往往无法通过直接的粒子级模拟来实现。

方法论
作者引入了一种基于约束的过阻尼布朗动力学方法，该方法用几何投影步骤取代了对刚性近接触力的显式时间积分。该方法作为 LAMMPS 模拟软件包中的一个 fix 实现，其概念灵感源自计算机图形学中使用的基于位置的动力学（position-based dynamics）。

该算法分为两个交替阶段运行：

1. 自由布朗更新：粒子在平滑的长程力（例如外部场、流体动力漂移）作用下进行标准的自由布朗动力学步骤，产生一个无约束构型 $X^{(0)}$。
2. 几何投影：将无约束构型投影到由几何约束定义的容许构型集上。这些约束强制执行：
   • 硬球排斥：非键合粒子不能重叠（$r_{ij} \geq d^0_{ij}$）。
   • 粘性键合：当粒子靠近至阈值距离内（$r_{ij} < d^0_{ij} + \epsilon$）时，形成不可逆键合并保持固定。
   • 硬壁相互作用：粒子受边界限制，保持最小间隙。

投影过程使用一种弛豫方案进行迭代求解，其中构型 $X^{(k)}$ 根据约束残差 $C_\alpha(X)$ 和迁移率矩阵 $M$ 进行更新。计算修正项 $\Delta X^{(k)}$ 以最小化残差，并采用欠松弛因子 $\omega < 1$ 以确保在稠密、键合网络中的收敛鲁棒性。键合形成通过更新约束集来处理，以确保投影容差与键合创建准则之间的一致性。

核心贡献

• 算法框架：一种用于过阻尼布朗动力学的模块化约束投影框架，将平滑相互作用（驱动自由步骤）与奇异硬/粘性相互作用（定义容许几何结构）分离。
• 实现：作为 LAMMPS fix 的工作实现，可通过 Zenodo 存档获取，允许模拟高密度受约束胶体系。
• 验证：通过针对两种基本极限的显式势能布朗动力学模拟验证了该方法：
  1. 非聚合硬球：投影动力学在体积分数 $\phi=0.55$ 和 $\phi=0.25$ 时重现了平衡径向分布函数和长时间扩散系数。
  2. 不可逆粘性聚合：该方法恢复了粘性硬球的聚合动力学，与具有日益短程吸引势的显式势能模拟结果相匹配。

结果

• 计算效率：基于约束的方法允许的时间步长（$\Delta t$）比显式势能布朗动力学大三个数量级（50 ns 对比 25–50 ps）。尽管投影求解器本身具有计算成本，但在测试的体积分数范围内，净物理时间吞吐量（$t_{phys}/t_{wall}$）增加了一到两个数量级。
• 电泳沉积 (EPD)：该方法被应用于模拟低 Peclet 数下的电泳沉积，这一领域此前因所需的漂移时间过长，在粒子级模拟中难以实现。
  • 在大约 1.5 天和 6 天内，分别实现了 $Pe$ 为 10 和 1 的模拟；而显式势能模拟则需要数周甚至无法实现。
  • 新的动力学见解：在高 $Pe$ 下，模拟证实了快速致密化以及随后由前沿诱导的键合过程。在低 $Pe$（$Pe=1$）时，观察到了截然不同的路径：悬浮液首先致密化，随后才发展出键合网络。这种低 $Pe$ 机制也导致最终沉积结构中形成的孔隙比高 $Pe$ 下发现的微小间隙孔隙更大。

意义
本文声称，当刚性近接触相互作用的物理效应可以用几何约束和接触捕获规则表示时，可以有效地从时间积分中移除这些相互作用。所提出的基于约束的布朗动力学方案保留了硬球结构和接触限制聚合的目标观测值，同时显著延长了可及的物理时间尺度。这种加速使得粒子级模拟低 Peclet 数下的电泳沉积变得可行，揭示了此前无法触及的沉积路径和结构特征。作者指出，该公式是模块化的，理论上可以扩展到非球形颗粒，但向可逆接触或依赖于构型的流体动力迁移率的扩展仍需进一步开发。主要成就在于将布朗动力学模拟的数值时间尺度与感兴趣的集体时间尺度对齐。