最近帮一个师弟调节分子动力学模拟的参数项目正好是矿物表面润湿性具体场景是把水滴放到石英表面上用模拟去算接触角。这个方向在浮选、油藏开发、土壤污染治理里都很重要实验上测接触角看着简单但界面附近的分子排列、水化层结构、前驱膜这些东西实验手段很难直接看到分子动力学模拟刚好能把这些细节拉出来。这篇文章我会把整个研究思路、建模流程、接触角计算方法还有一堆容易踩的坑都整理出来给正在做或者准备做矿物表面润湿性模拟的朋友参考。1. 为什么分子动力学模拟适合矿物表面润湿性研究1.1 实验看不全的界面结构MD能补上接触角听起来就是“液滴在固体表面形成的角度”实验上用测角仪测量几秒钟就能出结果。但有一个问题一直困扰着实验研究者液滴底部和固体表面之间存在一层只有几纳米厚的区域实验光学手段很难分辨这一层里水分子的密度和取向。比如石英表面到底是亲水还是疏水宏观上取决于表面的羟基浓度和氢键网络但实验很难直接“看见”氢键网络是怎么铺开的。分子动力学模拟恰恰是从分子坐标出发把每个水分子、每个表面原子的位置和速度都记录下来。你可以在模拟轨迹里直接看第一水层在表面上的密度可以看到水分子朝向表面排列还是背离表面排列还能统计氢键的数量和寿命。这些都是实验无法直接给出的微观机制。所以这个方法的项目定位很清晰宏观测角实验给结果分子动力学模拟给解释。两者结合才能回答“为什么这种矿物表面亲水那种表面疏水”这类问题。浮选药剂怎么吸附、油藏注水把这些机理搞清楚优化策略才做得动。1.2 润湿性模拟里的核心观测量在分子动力学模拟里我们通常关注的观测量有四个接触角。这是最直接的润湿性指标通常从液滴密度轮廓中提取。固液界面张力。通过杨氏方程把接触角和界面能联系起来但直接计算界面张力比较麻烦。水分子密度分布。特别是垂直表面方向的密度剖面能看出是否存在致密水化层。水分子取向和氢键网络。能解释表面亲疏水的分子机制。接触角的定义来自杨氏方程cosθ (γ_sv - γ_sl) / γ_lv在模拟里要算接触角有两条路线。一条是直接模拟液滴让液滴平衡后从密度场里读出轮廓再拟合得到角度另一条是分别算三个界面张力代入杨氏方程。第一条路线简单直接也是大多数人采用的办法第二条路线对力场和静电处理的要求高如果体系带大量电荷界面张力收敛很慢容易算偏。你不需要一开始就把所有量都算出来。做矿物润湿性模拟先把接触角算稳再补充密度剖面和氢键统计基本就能构成一篇工作的核心数据了。2. 建模与力场搭配矿物表面润湿性模拟的第一步2.1 矿物表面模型怎么搭才靠谱我以最常见的石英为例。石英的晶胞参数已知沿着特定晶面切一刀比如石英(001)面然后用晶格参数扩展成超胞表面尺寸大约 6×6 纳米。这样的大小对于几纳米的水滴来说基本够用周期性镜像的影响也可以接受。切表面以后要处理悬键。裸露的硅原子需要加上羟基-OH模拟中一般把这层羟基视为接近真实石英表面水化层的状态。有人会问不加羟基行不行实验结果已经表明干燥石英表面的润湿性和水化后的石英完全不同。如果不做羟基化接触角会被严重高估算出的表面可能变成疏水这和实际情况完全相反。具体操作上你可以用 Materials Studio、GULP 或者 VMD 里的无机结构构建脚本。我的习惯是在 Materials Studio 中完成表面切割和羟基化再导出为 PDB / data 文件导入到 LAMMPS 中继续做。对于黏土矿物常见做法是把层状结构切开在层表面补上羟基和层间阳离子再调整晶胞尺寸让表面接近中性。有一点必须注意石英表面的质子化状态取决于 pH。酸性环境下表面羟基以 Si-OH 为主中性到碱性环境下部分会解离成 Si-O⁻。如果你模拟的是中性水环境设定成全羟基化是可以接受的但如果你要研究 pH 对润湿性的影响就需要显式改变表面电荷和抗衡离子不能简单只用中性表面。2.2 水模型与力场的匹配问题矿物表面润湿性模拟里的力场选择基本决定了结果的可靠性。水分子模型常用 SPC/E、TIP3P、TIP4P三种模型对体相水密度、扩散系数、表面张力的描述略有差异。对于石英和黏土矿物我建议优先考虑 CLAYFF 力场因为它本身就是针对矿物与水界面开发的氧、硅、铝、镁等参数和 SPC/E 水兼容性很好。常见的组合是 CLAYFF SPC/E。这个组合对石英表面接触角模拟的复现性不错是我目前用得最多的一套。另外一个选择是 Interface 力场Interface FF它对有机分子修饰矿物表面的体系更合适尤其是表面存在有机官能团的时候。力场文件中的非键参数特别是短程排斥参数和原子电荷直接决定水分子在表面上的吸附强度。如果参数太强第一水层会过于致密接触角偏低如果参数太弱水分子吸附不上去接触角又变高。所以拿到一个力场组合要先算一个验证体系把模拟得到的体相水密度、径向分布函数和实验值对比再算石英表面的水化层结构和实验X射线反射数据对比。验证通过后再去做液滴接触角才不会浪费机时。2.3 力场参数选错的典型后果我做过的几次失败经验里印象最深的是选了一组不匹配的力场表面羟基沿用 CLAYFF但水分子用了 TIP4P。看似一个细节实际接触角比 SPC/E 组合偏高了 15°左右。原因是 TIP4P 对水的表面张力描述和 CLAYFF 的非键参数在界面混合规则上存在偏差导致水在表面的铺展能力下降。还有一次是在模拟方解石表面时忽略了方解石表面钙离子的电荷平衡。方解石表面原本带正电需要配一定的碳酸根或者抗衡阴离子但我直接用了电中性条件结果水分子在表面上的排列毫无规律接触角波动巨大。后来补上表面碳酸根化学计量比正确以后结果才稳定下来。所以给一个硬性建议在模拟化学计量比和电荷平衡出现问题之前先用一个 3 纳米见方的小体系做 200 皮秒的测试模拟观察水分子的密度剖面是否合理。如果发现水分子远离表面或者紧贴表面成为一层冰状结构就要回头检查力场参数和表面电荷。3. 实操流程从初始构型到接触角计算3.1 初始液滴摆放和模拟参数一个比较稳妥的模拟流程是这样生成石英表面超胞厚度至少 1.5 纳米xy 方向尺寸 6 纳米以上。在表面上放置一个预先平衡好的球形水滴半径设 3 纳米左右大约包含 3000~4000 个水分子。水滴初始构型可以直接在 VMD 里生成一个水盒子再切成球形也可以用 packmol 生成球形区域内的水分子。把水滴放到距离表面 0.25 纳米左右的位置。不要直接贴在表面否则初始势能太大水分子会“弹开”。设置周期性边界条件z 方向真空层高度要比水滴直径多 2 纳米以上防止水滴和周期性镜像相互作用。选择 NVT 系综温度设置为 300K用 Nosé-Hoover 恒温器时间步长 1 飞秒。先用能量最小化消除局部接触再跑 1 纳秒平衡让水滴在表面上铺展到稳定状态随后采集 3 到 5 纳秒的轨迹用于分析。这里有个细节如果初始水滴离表面太远需要很长时间吸附如果太近水分子的氧会和表面硅原子产生强烈排斥直接导致模拟崩溃。我一般用软势或者先跑几步 NVE 让体系自己调整然后再切到 NVT 稳定温度。另外表面原子最好固定不动。虽然表面原子弛豫在真实体系里也存在但在润湿性模拟中如果表面原子参与热运动接触角结果会变得很不稳定。通过 fix spring 或 velocity zero 处理表面原子既能保持表面结构又能避免表面层整体漂移。3.2 接触角提取的密度剖面法从模拟轨迹里提取接触角最常用的是密度剖面拟合法。思路很简单把圆柱坐标以水滴质心为轴在 xy 平面内划分同心圆环每个环的宽度约 0.2 到 0.5 纳米。对每个圆环统计沿 z 方向的水密度分布得到水密度随径向距离 r 和高度 z 的变化。找到液滴轮廓的等密度面一般取体相水密度一半作为分界线。体相水密度在 SPC/E 模型下约为 1 克每立方厘米对应数密度约 0.033 个每立方埃。把等密度面上的点提取出来拟合球冠方程得到球冠的半径和中心高度进而计算接触角。拟合球冠时可以用如下近似关系如果液滴的球冠半径为 R冠状高度为 h底面半径为 r0那么接触角 θ 满足cosθ (R - h) / R如果 h R则 θ 小于 90°表面偏亲水如果 h 接近 R则 θ 接近 180°表面疏水。实际操作中等密度面靠近固体的部分会受到吸附水层干扰。石英表面通常有一层密度很高的水化层这部分水属于表面吸附水不算液滴主体。如果你不把这一层排除拟合出来的接触角会偏小。我的做法是先算垂直表面的密度剖面找到第一水层的密度峰位置然后在等密度面提取时忽略表面 0.5 纳米以内的区域只拟合液相主体轮廓。3.3 影响接触角数值的几个隐藏因素接触角并不是一个完全稳定的“物理常数”在分子动力学模拟里受几个因素影响很大。第一个是温度。温度升高会增加热涨落液滴形状波动变大拟合误差自然增大。更头疼的是如果真空层里面水分子饱和蒸气压比较高水分子会蒸发到真空层里形成少量气相水虽然对密度等值线影响不大但会干扰质心判断。我一般把温度控制在 280K 到 320K 之间并且用 NVT 而不是 NVE让温度稳定在目标值。第二个是液滴尺寸。纳米尺度下的接触角存在明显的线性尺寸效应类似杨氏方程的修正形式cosθ_app cosθ∞ - τ / (γ_lv · R)其中 τ 是线张力R 是液滴底半径。液滴越小该项修正越大。因此模拟里不能只算一个尺寸就下结论至少要分别模拟半径 2.5、3、4 纳米的水滴把接触角随尺寸变化趋势画出来外推到无限大液滴得到宏观接触角。第三个是采样时间。很多人跑 500 皮秒就去算接触角结果误差很大。液滴在表面上铺展是一个慢过程尤其是水滴在疏水表面时铺展速度更慢。我建议至少平衡 1 纳秒采样 3 纳秒以上并且把时间窗口分成两段分别计算接触角确认两段结果一致再取平均。4. 常见问题与排查接触角飘、水滴飞、静电截断4.1 水滴在表面“跑”了怎么办模拟中经常出现水滴逐渐偏离表面中央甚至跑到盒子边缘的情况。表面如果完全均匀这种情况不太明显但如果你用了带缺陷的表面或者无定型表面局部势能不均匀水滴就会向势能更低的地方迁移。解决方法有几个固定表面原子后水滴的质心漂移大部分来自初始摆放位置不在中心。放置液滴前把表面质心和水滴质心都投影到 xy 平面并调整到同一位置。可以在 xy 方向加上一个非常弱的约束势比如对水滴质心加一个弹性常数很低的固定不影响接触角但能防止水滴跑出采样区域。检查周期性镜像。如果盒子 xy 尺寸只有液滴尺寸的 2 倍水滴的镜像会相互吸引导致液滴被“拉”向边界。我试过最省事的办法是在 LAMMPS 里给表面原子设置 position 固定然后每隔 500 步输出一次水滴质心坐标及时监控漂移。如果漂移速度超过每纳秒 0.5 纳米多半是表面电势分布有异常。4.2 接触角随时间一直漂移如果你把接触角按时间段分段统计发现数值不断下降或上升说明体系还没平衡或者表面结构在慢慢重构。无定型石英表面或者有表面缺陷的情况下表面原子可能缓慢弛豫导致表面能变化接触角跟着漂移。解决办法是先单独对表面进行 1 纳秒以上的松弛确认表面构型稳定后再放水滴。另外如果表面羟基发生质子转移虽然 CLAYFF 中不会自动发生但某些反应力场例如 ReaxFF 会出现接触角也会漂移。还有一种情况是水滴用的太少例如只有几百个水分子热涨落引起的接触角波动可能超过 10°。这时候增加水分子数量比延长模拟时间更有效。4.3 盒子尺寸和周期性镜像效应周期性盒子虽然方便计算但会带来两个麻烦一是水滴在 xy 方向上的镜像会互相影响二是 z 方向真空层太薄会导致水滴和上方镜像“隔空作用”。盒子尺寸的经验值xy 尺寸至少是液滴底直径的 3 倍。例如半径 3 纳米的水滴底面直径约为 5 纳米盒子 xy 方向至少做 15 纳米其实有点夸张了。实际操作中 6~8 纳米也可以但必须做尺寸收敛测试把盒子边长从 6 纳米增加到 10 纳米看接触角变化是否超过 2°。z 方向的真空层高度建议至少为液滴直径的 1.5 倍。因为静电是长程作用即使有 PPPM周期性镜像之间的静电耦合也不能完全消除。LAMMPS 里处理这个问题的标准做法是设置kspace_modify slab 3.1092这样会给 z 方向一个校正项近似屏蔽周期性镜像。如果不用这个 slab 修正你会看到水滴在 z 方向上出现一层不自然的高密度薄片。4.4 长程静电处理不当接触角直接不可信矿物表面和水都是带电体系如果不做长程静电处理结果会非常离谱。简单 cutoff 截断会让水分子的偶极-偶极相互作用产生截断伪像导致表面附近的水结构偏离真实状态。我用 LAMMPS 时的标准配置如下pair_style hybrid/overlay lj/cut/coul/long 10.0 pair_coeff * * lj/cut/coul/long 0.0 0.0 kspace_style pppm 1e-4其中 LJ 截断设为 10 埃静电用 PPPM精度 1e-4。如果你用的是 CLAYFF这是比较稳妥的配置。一个容易犯的错误是把真空层选得很大然后三维周期性 PPPM 会在 z 方向产生一个不需要的电场。要解决这个问题除了加 slab 修正也可以把 z 方向设为非周期性但那样就不能用标准的 PPPM 了。我更推荐保留周期性晶体 slab slab 修正这也是多数文献里用的做法。5. 扩展思路从简单水润湿到油水矿物三相体系5.1 把水滴换成油水双相做矿物表面润湿性的研究很多时候不只是水-空气-矿物三相而是油-水-矿物三相。尤其在工作机制和浮选机理分析里我们要知道矿物表面在油藏条件下到底亲油还是亲水。这时可以把体系扩展成矿物表面上方先铺一层水再铺一层有机液体比如正癸烷然后观察油水界面在矿物表面上的形态。这个扩展看着不难但有几个坑油分子和矿物的力场参数需要单独验证。正构烷烃可以用 TraPPE 或 OPLS-AA但跨力场混合规则是否和 CLAYFF 兼容需要先算油在矿物表面的吸附能对比实验值。油水界面张力也很容易算偏。需要在模拟前先算一个油水界面体系验证界面张力在合理范围内。三相接触角的提取变得更复杂。不能只用水密度轮廓还需要同时跟踪油分子和水分子的位置分布然后根据界面交线来确定接触角。我实际跑过油水石英体系发现水分子在石英表面形成一层非常稳定的水膜油分子很难直接吸附到表面。这和实验上石英表面亲水、油膜不稳定是一致的。如果模拟结果出现油直接贴附表面多半是力场混合参数高估了油和表面羟基的相互作用。5.2 我的几点实操体会做矿物表面润湿性模拟这几年我最大的体会是先别急着追求复杂的模型。先用一个简单的石英表面加水滴把接触角算到和实验值一致再往上加有机分子、修饰表面、改变 pH这样一步步推进最容易排查问题。另外轨迹存储频率别太低。接触角分析需要足够多的液滴轮廓帧数建议每 1000 步保存一帧这样可以计算误差棒。如果只保存最后几帧误差分析就无从谈起。还有一个小技巧在跑正式模拟之前先写一个分析脚本用一个人造液滴轮廓例如已知接触角为 60° 的半球冠测试脚本输出的接触角是否正确。这个步骤能避免你花了两周机时最后发现是脚本代码写错了。最后如果你希望发表高质量的结果建议除了接触角至少补充一个界面微观证据比如水分子在表面附近的密度峰位置或氢键数量。这样实验上看到宏观接触角差异模拟里能找到分子层面的原因整个工作才算闭环。