- 原文 Blowing Up
- 2016-12-22 15:16:04; 翻译: 刘恒江; 校对: 李继存
【刘恒江 按】我们使用GROMACS进行分子动力学模拟时常常会遇到LINCS warning这样的警告, 过多的警告会导致体系崩溃, 程序运行异常. 出现LINCS warning往往意味着体系初始构型不够合理, 在模拟过程中出现了异常键. 对于这种情况的应对方法, GROMACS官网上其实已有了详细的介绍, 在这里翻译一下, 供大家参考.
爆破(Blowing Up)是一个用于表述模拟失败的相当专业的术语. 简要地说, 它描述了一个由于产生了极大的作用力并最终导致积分运算失败的典型错误.
稍微展开一下背景, 我们必须清楚分子动力学的基本原理是在非常短的, 离散的时间步长内对牛顿运动方程进行数值积分, 并借助这个时间步从粒子前一时间步的速度, 位置, 受力来确定下一时间步的速度和位置. 如果在一个时间步中作用力变得很大, 这将导致粒子在到达下一时间步时会在速度/位置上产生极大的变化. 典型情况下, 这将激发一连串致命的错误: 一个原子在某一时间步中受到了很大的作用力, 在下一步中它可能失控并击穿整个体系, 最终停在超出它应处的范围, 与其他原子重合的位置, 或类似的情形. 这又使得下一时间步中产生了更大的作用力, 更不受控的运动. 这种情况延续下去, 最终, 这将导致模拟程序在某些方面崩溃, 因为它没法处理这种情况. 在使用了约束的模拟中, 这个问题最初的征兆通常就是出现一些LINCS或SHAKE的警告或错误——这并不是因为这些约束导致了错误, 而是因为它们最先受到影响而崩溃. 类似的, 在模拟中使用区域分解时, 你可能会看到类似粒子处于其电荷组区域分解区外超过一个区长度的信息, 这也是你的体系潜在问题的征兆, 而不是DD算法自身的问题. 此外, 有关表格或1-4相互作用处于表格支持范围外的一些警告也是一样的. 由于这些计算在不同计算机系统间不具有数值重现性, 一台计算机中出现问题的模拟在另外的计算机中可能会获得稳定的体系.
导致爆破问题的可能原因有:
- 能量最小化不彻底;
- 初始构型不合理, 可能有空间冲突;
- 采用了过大的时间步长(尤其是在使用约束的时候);
- 对自由能计算中的粒子插入没有使用软核;
- 使用了不正确的压力耦合算法(例如: 在还没有平衡时, Berendsen方法能够很好地弛豫体积, 但接下来就需要换到更精确的压力耦合算法了);
- 使用了不正确的温度耦合, 或许用在了不正确的组上;
- 对粒子坐标进行的位置限制与当下体系中的坐标差别太大;
- 体系内某处有一水分子与其他水分子独立开来;
- 遇到了
mdrun中的bug;
由于爆破是由于特定时间步长下作用力过大导致的, 最基本的解决方法有以下两种:
- 确保作用力不会过大;
- 使用更小的时间步长;
如果问题出现在模拟的开始阶段, 更好的体系准备工作将有助于第一条.
如何诊断一个不稳定的体系
对一个爆破的体系进行错误排查是一项极具挑战性的工作, 尤其是对那些分子模拟的萌新们来说. 当遇到这种情况时, 以下几点建议通常会有所帮助:
- 如果程序崩溃发生得相当早(几步之内), 设置
nstxout(或nstxtcout)为1, 收集所有可能的帧. 观察结果轨迹去排查是哪个原子/残基/分子首先变得不稳定; -
逐渐简化体系来寻找原因:
- 如果你模拟的是一个溶剂构成的盒子, 试着在加入溶剂前先对单个溶剂分子进行能量最小化以及模拟, 查看体系的不稳定性是否是由分子拓扑的内在问题引起的, 或者在初始构型中分子位置是否有冲突;
- 如果你模拟的是蛋白质-配体体系, 试着将蛋白质单独放置在所需溶剂中进行模拟. 如果蛋白质是稳定的, 将配体分子单独放在真空中进行模拟, 观察其拓扑是否能提供稳定的构型, 能量等;
- 删除设定的算法(LINCS/SHAKE), 尤其在没有达到充分平衡的时候;
- 使用
g_energy(5.0版本后为gmx energy)来监测体系各个能量组分的变化, 如果分子内的能量变化出现了尖峰, 那可能意味着使用了不正确的键参数; - 确保你没有忽略在运行
mdrun之前的步骤中的任何错误信息(运行pdb2gmx时缺失原子, 运行grompp时原子名称不匹配, 等等)或者使用了变通方法来保证你的拓扑文件完整并能被正确地解释(如在grompp时使用了-maxwarn选项忽略不能忽略的警告); - 确保在
.mdp文件中对你的体系和所选力场设置了正确的选项, 尤其是截断的处理, 近邻搜索间隔(nstlist)以及温度耦合. 就算体系的初始构型是合理的, 不正确的设置也会导致物理模型的崩溃;
如果使用隐式溶剂, 在预平衡阶段使用比成品MD更小的时间步长可以使能量分配更稳定.
在一些常见的情形下不稳定现象出现的频率很高, 通常是向体系中引入新物种(配体或其他分子)的时候. 想要找到问题的根源, 就得将体系拆开来逐步分析(如蛋白质-配体复合物体系):
- 蛋白质自身(在水中)能否充分进行能量最小化? 这是对蛋白质分子坐标完整性和体系准备工作的一个测试. 如果失败了, 说明在运行
pdb2gmx中可能出错了(说明见下面), 或者使用genion添加离子时, 离子的位置离蛋白质太近了(毕竟它是随机加入的); - 配体分子在真空中能否进行能量最小化? 这是对其拓扑的测试. 如果不能, 检查你对配体分子的参数化设置以及任何加入力场文件中的新参数;
- (如果第二步中的能量最小化成功)配体分子在水中能否能量最小化, 或直接运行一段短时间的模拟?
还有一些问题可能来自于生物分子的拓扑自身:
- 在运行
pdb2gmx时是否使用了-missing选项? 如果是, 请删除该选项. 重新构建缺少的坐标而不是忽略它们; - 是否改变键长并不顾长/短键的警告? 如果是, 请不要这么做. 这将导致原子缺失或形成一些糟糕的输入结构.
【刘恒江 补充说明】
总得来说, LINCS/SHAKE warning的出现还是因为体系没有平衡好. 往往我们在进行了能量最小化后还是偶尔会遇到这个问题, 而且一般出现在NPT平衡阶段(NVT因为没有加压所以出现警告的概率较小). 在确认拓扑没有问题的前提下, 可以先进行一段时间步长较小的平衡来松弛体系, 往后再进行较大时间步长的平衡. 当然, 也可以使用不同的压力耦合算法来实现快速松弛体系的效果.
在GROMACS邮件列表里面遇到这个LINCS warning的人非常多, 同样列出些大家分别遇到的问题, 以便参考:
