刘建川・GROMACS 实战 | 计算液体粘度：cos-acceleration怎么用？
大家好，我是刘建川，14 年深耕 GROMACS 分子模拟领域，始终专注分子模拟实战。在分子动力学模拟中，密度、RMSD、扩散系数这些指标大家接触得比较多，但如果研究的是水、有机溶剂、离子液体、电解液等液态体系，粘度（Viscosity）同样是一个很重要的性质。
一方面，它可以用来和实验值进行比较；另一方面，粘度对分子间相互作用比较敏感，也能从侧面反映力场对体系动力学性质的描述是否合理。
GROMACS中计算粘度的方法不止一种。本文主要介绍一种比较经典的做法：
周期性微扰法，也就是 cos-acceleration 方法。
一、基本原理
这个方法的思路并不复杂。
在模拟过程中，沿 x方向 给体系施加一个随 z坐标 周期变化的加速度：
其中：
A：施加的加速度 Lz：模拟盒子在 z 方向的长度
经过一段时间后，液体会形成相应的速度分布：
最终可以根据加速度、速度响应、体系密度以及盒子尺寸计算粘度：
实际使用时并不需要自己重新实现这些计算，GROMACS已经提供了相应功能。
二、模拟体系怎么准备？
建议先把体系正常平衡好，再进行粘度计算。
常见流程可以写成：
体系构建 → EM → NVT → NPT → NVT → 粘度模拟
前面的 EM、NVT 和 NPT 与普通分子动力学模拟没有太大区别。
这里主要注意两个问题。
第一，密度要基本稳定。
因为粘度计算本身涉及体系密度，如果 NPT 阶段体系还没有平衡好，后面的结果自然也不太可靠。
第二，粘度模拟阶段建议使用固定盒子尺寸。
所以通常可以在 NPT 平衡结束后，选取稳定结构，再进行一段 NVT，然后进入后续非平衡模拟。
此外，这种方法中经常会把 z方向适当拉长，例如：
4 × 4 × 12 nm³
原因是速度场沿 z 方向变化，适当增加 z 方向长度通常有利于获得更清晰的速度响应。
注意：盒子尺寸并没有一个固定比例，具体还需要根据体系大小和计算量来调整。
三、关键参数：cos-acceleration
真正进行粘度模拟时，在 .mdp 文件中加入：
cos-acceleration = 0.025
单位为：
nm ps^-2
这个参数决定施加给体系的周期性外部加速度。
这里有一个比较容易踩的坑：
0.025并不是固定值，也不是所有体系都适用。
如果加速度太小，产生的速度信号可能被热涨落淹没；如果加速度太大，又可能让体系偏离正常的线性响应范围。
实际计算时，可以尝试几个不同的值，例如：
cos-acceleration = 0.01
cos-acceleration = 0.025
cos-acceleration = 0.05
然后比较不同条件下得到的粘度。
如果随着外场增大，计算得到的粘度明显下降，就需要考虑剪切强度是否已经过高。
四、模拟结束后怎么看结果？
模拟结束后，可以直接查看能量文件：
gmx energy -f vis.edr
开启 cos-acceleration 后，GROMACS会记录体系的速度响应，同时输出和粘度计算有关的数据。
从操作上来说，这一步其实并不复杂。
真正需要注意的是：算出一个粘度值，并不代表这个结果一定可靠。
五、粘度结果至少检查这3件事
1. 模拟时间是否足够
粘度属于波动比较明显的动力学性质。
如果模拟时间太短，很容易出现前几纳秒一个结果、后几纳秒又发生明显变化。
因此最好比较一下不同时间段的平均值，看结果是否逐渐稳定。
2. 不同加速度下结果是否一致
cos-acceleration 不能设置得过大。
如果条件允许，建议至少测试 2～3个不同的加速度。
如果在较小外场范围内得到的粘度比较接近，结果通常会更有说服力。
3. 不要忽略力场本身
这一点很重要。
有些力场可以很好地描述液体密度，却不一定能够同时准确描述扩散系数和粘度。
因此，如果模拟粘度和实验值存在一定偏差，不一定是模拟流程出了问题，也可能与力场参数本身有关。
六、简单总结
使用 GROMACS 计算液体粘度，大致可以归纳为：
先充分平衡体系，再利用 cos-acceleration 施加周期性外场，通过液体产生的速度响应计算粘度。
实际操作并不算复杂，但有三个问题值得重点关注：
体系有没有平衡好、外场是不是太强、结果有没有真正收敛。
相比最后得到一个具体数字，这几个问题往往更重要。
如果你正在做水、有机溶剂、离子液体、电解液或其他液态体系的 GROMACS 模拟，这个方法可以作为粘度计算的一种参考。
