Martini力场在体系自组装时报错-nan和too many LINCS warning

zerochan

（系统为ubuntu22.04, Gromacs2022.02）最近使用Martini力场（2）进行多肽自组装（多肽可能较大，有上百个氨基酸），按照官网上多肽自组装教程，在进行加水替换离子之后对体系进行最小化后，发现体系在设置的5000步条件下，差不多500多步就结束了，提示力无法收敛到最小值或以达到机器精度（如图1），然后进行平衡的时候开始没几步就一直提示LINCS WARNING，后面步骤就有些没有这个LINCS WARNING，很多有LINCS WARNING，最后体系就崩溃了，报错Fatal error:Too many LINCS warnings (1103)等等（见图二）。根据网站上一些信息和手册，调整了平衡mdp文件中的dt，nsteps，lincs_order，lincs_warnangle来调整积分条件和约束限制。但是随后出现了警告 Warning: pressure scaling more than 1%, mu: inf inf inf  

Warning: Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.

Box (3x3):
            Box[    0]={         inf,         -nan,         -nan}
            Box[    1]={        -nan,          inf,         -nan}
            Box[    2]={        -nan,         -nan,          inf}
         Can not fix pbc.
Internal error (bug):
Step 100: The total potential energy is -nan, which is not finite. The LJ and
electrostatic contributions to the energy are 0 and 0, respectively. A
non-finite potential energy can be caused by overlapping interactions in
bonded interactions or very large or Nan coordinate values. Usually this is
caused by a badly- or non-equilibrated initial configuration, incorrect
interactions or parameters in the topology.
然后体系又崩溃了（见图三）。

根据提示，可能是体系不平衡，拓扑结构有问题，随之又修改了最小化的mdp文件的dt，nsteps，lincs_order，lincs_warnangle等等参数，还进行了steep和cg等多步最小化，但是力一直收敛不了，最后得到的最小化文件跑平衡，会出现上述的问题，有时显示rmsd特别大，有时到系统崩溃rmsd还是在0.00几，但是一直解决不了。所以准备检查文件结构是不是有问题，因为水是直接从网上获得的，所以没检查，在检查单个蛋白结构的时候，开始跑蛋白在水中结构的时候，前几步也出现了 LINCS WARNING，但是到十几步的时候就停住了，体系没有崩溃，但是跑了好几天了一直都停在那里，也不知道是卡住了还是在继续跑（见图四）。

所以出现如下问题：
1.对于上述问题，是不是拓扑的问题？我该如何检查和解决？现在没崩溃的体系还在跑吗？要停止还是继续等？

2.对于mdp里头单个参数粗略知道含义，但是它们之间的联系和限制有很多不清楚，看手册只知道rlist要>=rcoulomb和rvdw，但对其他的还是搞不清它们之间的联系，请问是手册那一章说明这些或者有什么资料可以进行学习的？

3.在Martini自组装教程中，他们可以算出这个盒子中要添加多少水分子，这个体系中肽的浓度，这些是根据什么条件确定的呢，这个计算是不是只是根据盒子大小，加入肽分子的量，用普通的摩尔计量算出来的（因为我在上述体系加水过程中，会发现有些时候这个水在盒子中会有叠加，盒子中水密度不同）？

4.对于化学计算的很多参数，特别是时间，假如我实验条件用的10mM，30min，不同的加热温度或物理场，这些在计算机模拟下如何体现，时间这么长实现肯定不现实，有些物理场比如等离子也模拟不了，这样如何对等实验和计算？是不是只能说现在计算只是一个参考，无法真正的指导实验？

哈哈哈，因为非量化专业，最近也在恶补这些知识，理论实践水平有点低，只能根各位专业的同学老师请教一下，分场感谢！！！

附上后期最小化和平衡的mdp，再次谢谢各位！！！

Frozen-Penguin

zerochan 发表于 2022-10-23 10:11
非常感谢，其实现在试了一下，包括dt（1fs），tau_p, tau_t和能量最小化（生成蛋白CG，水，多个蛋白，多 ...

能量最小化的Fmax偏大也能继续模拟，而且体系崩溃的时候已经模拟很长时间了，基本可以排除能量最小化没做好的可能，引起误差的原因可能是控压，用gmx energy查看模拟过程中的压强变化，如果是在崩溃前一小段变化幅度突然变大，就应该继续改大tau_p和改小dt。
软核势可以用在计算自由能中，所以手册把两者写在一起了，lambda和delta G都是自由能相关的内容，软核势可以单独用，看手册中的示意图设置参数就可以了。

zerochan

Frozen-Penguin 发表于 2022-10-14 23:24
1.能量最小化的时候最大的力有些偏大了，可能是力场不合理，需要检查，可以尝试加Maritni的弹性网络，增 ...

非常感谢，我的具体问题在回复q11pet28abl21朋友时描述了，还是有问题，不知您是否清楚？还是非常谢谢，回头真的要好好学习下理论才行。

zerochan

q11pet28abl21 发表于 2022-10-15 20:40
之前在做肽段自组装遇到过这个情况。
我是使用steep再用CG去做能量最小化，再进行模拟时候出现这样的报错 ...

非常感谢，其实现在试了一下，包括dt（1fs），tau_p, tau_t和能量最小化（生成蛋白CG，水，多个蛋白，多个蛋白加水体系都应用了最小化），虽然可以跑短时间的平衡,但是长时间平衡还是会报错如下：
Step 1889000  Warning: pressure scaling more than 1%, mu: inf inf inf
step 1889000, will finish Sat Oct 22 09:39:53 2022Warning: Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.
         Box (3x3):
            Box[    0]={         inf,         -nan,         -nan}
            Box[    1]={        -nan,          inf,         -nan}
            Box[    2]={        -nan,         -nan,          inf}
         Can not fix pbc.
.....
Warning: Only triclinic boxes with the first vector parallel to the x-axis and the second vector in the xy-plane are supported.
         Box (3x3):
            Box[    0]={        -nan,         -nan,         -nan}
            Box[    1]={        -nan,         -nan,         -nan}
            Box[    2]={        -nan,         -nan,         -nan}
         Can not fix pbc.


-------------------------------------------------------
Program:     gmx mdrun, version 2022.2
Source file: src/gromacs/mdlib/sim_util.cpp (line 554)
Function:    void checkPotentialEnergyValidity(int64_t, const gmx_enerdata_t&, const t_inputrec&)

Internal error (bug):
Step 1889100: The total potential energy is -nan, which is not finite. The LJ
and electrostatic contributions to the energy are 0 and 0, respectively. A
non-finite potential energy can be caused by overlapping interactions in
bonded interactions or very large or Nan coordinate values. Usually this is
caused by a badly- or non-equilibrated initial configuration, incorrect
interactions or parameters in the topology.

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
看了一下错误信息和手册对应的问题，提示也是改上面的参数，或者本身结构有问题，
现在对自己这个体系猜测可能有一下问题在压力耦合我用的是Pcoupl = C-rescale 的选用与tau_p=12有问题（因为先前跑的时候提示BBerendsen和parrinello-rahman都不稳定，最后选的C），
还有就是结构有问题，导致能量最小化一直收敛不到体系要求的Fmax=100,总有部分原子为10^3-5次幂，跟错误信息中就是相邻原子（如2663,2664）即分子内结构有关（可能范得华力问题），手册提醒要用软核势，看了一下属实没看懂（什么Lamda，delta H），
因此想请问，如何解决这些问题，特别是结构上的优化和检查（软核势的应用，网上关于这个点的连接也不多），使最后体系能够跑到微秒级别？一直卡在这个问题上，或者有一些推荐的教学书或视频什么的？

q11pet28abl21

之前在做肽段自组装遇到过这个情况。
我是使用steep再用CG去做能量最小化，再进行模拟时候出现这样的报错。
我使用的方法是调小步长去跑5-10ns，等体系致密了再调回步长去正式跑，这时就不会出现报错。

Frozen-Penguin

zerochan 发表于 2022-10-14 15:37
非常感谢Frozen-Penguin的回答，
1.我此前试过单链加了水后进行多次能量最小化过，还是会出现上面跑多个 ...

1.能量最小化的时候最大的力有些偏大了，可能是力场不合理，需要检查，可以尝试加Maritni的弹性网络，增加系统的稳定性，但是这样可能不适合做组装，因为加了弹性网络后结构不会发生较大的变化
3.可能只是显示问题，也可能是因为水分子太多导致编号超过预留的位数导致数据错位，一般用gromacs处理gro文件是不会这样的，因为编号超过上限会归零，但如果是这样，这就可能是导致体系崩溃的原因，可以检查最后几个水的数据是否有问题，如果没问题就不用管显示结果
4.盒子的体积是知道的，如果知道摩尔浓度，自然可以算出应该加多少肽链

zerochan

Frozen-Penguin 发表于 2022-10-10 14:51
1.可以先做单链的模拟，在加水和离子之前先做一次能量最小化以避免因为粗粒化过程不合理而出现问题，模拟时 ...

非常感谢Frozen-Penguin的回答，
1.我此前试过单链加了水后进行多次能量最小化过，还是会出现上面跑多个单链体系的同样问题，最近对一开始生成单链蛋白结构进行最小化，对水单独进行最小化，合并水和单链蛋白后再做最小化，确实会跑的更就，但是最后结果还是一直体系崩溃。
在所有在最小化过程中都会有类似提示
Energy minimization has stopped, but the forces have not converged to the
requested precision Fmax < 100 (which may not be possible for your system).
It stopped because the algorithm tried to make a new step whose size was too
small, or there was no change in the energy since last step. Either way, we
regard the minimization as converged to within the available machine
precision, given your starting configuration and EM parameters.

Double precision normally gives you higher accuracy, but this is often not
needed for preparing to run molecular dynamics.
You might need to increase your constraint accuracy, or turn
off constraints altogether (set constraints = none in mdp file)

writing lowest energy coordinates.

Steepest Descents converged to machine precision in 1969 steps,
but did not reach the requested Fmax < 100.
Potential Energy  = -1.6663042e+08
Maximum force     =  1.9552879e+04 on atom 1098
Norm of force     =  3.1101550e+01
指出力无法收敛到指定大小，但是符合机器精度，看了一下这个emtol可以设置，因为最小化的结果总会有原子和默认的相差两个以上数量级，这个力要是否要根据体系大小来设置，依据是什么？而在体系崩溃的时候起初前几步会有LINCS WARNING，键的旋转超过了预设值，最后体系崩溃提示Fatal error: Bond length not finite.
dt我都是用的2fs，压力温度tau我调的都很大了压力我调到了12,温度调到了3，最小化后用pymol查看发现这个大多肽或者说是蛋白有些都不成一个整体了
猜测是不是在每一步约束上要对H原子进行约束，约束没设置好？还是体系太大不加水单链有几千个原子，插入蛋白单链有几万个原子，加水后有上百万原子。我电脑是AMDR9，1650Mq，现在一般这种聚集最好是低于多少原子会比较稳妥？

2.嗯嗯，后续会继续深化手册学习，每个参数代表意义大概知道代表什么，不过里头涉及的联系，特别是数学上的还是学起来有点难，就是相关多个参数之间的联系，短期内感觉要把基本几百页的说明书看完感觉到看懂确实能力有限。只能继续加紧学了。

3.我上传一个新图，里头加水后，盒子的水在不同区域密度不同，不太了解为啥没有平均？

4.可能我没表达清楚，就是在教程中，他们多肽浓度可以计算出来，在这些定量关系下要如何定义相对大小关系？ 最后倒数第二图，等了一个星期还没结束，不知道是卡住了还是还在跑，很是迷惑？我该手动停了吗？


哈哈，非常感谢，确实不会的太多，非常感谢，现阶段就差这步了，唉

Frozen-Penguin

1.可以先做单链的模拟，在加水和离子之前先做一次能量最小化以避免因为粗粒化过程不合理而出现问题，模拟时改小dt或者改大温度压强平衡的时间常数tau，这样可以使体系相对更稳定。
2.学习mdp各参数的含义需要先了解分子动力学模拟具体是怎么实现的，大体了解之后看手册就能看懂，一般只需要知道几个重要的参数怎么设置就行了，大部分参数都可以用默认的。
3.加水的原理是复制一个经过平衡的水盒子，填满整个系统，然后把和体系中已有的分子重合的水分子删掉。
4.粗粒化模型的时间和温度一般都无法和现实的时间和温度严格对应，只能反映相对的大小关系，所以这些模拟只能给出一个趋势，没有必要追求准确的数值。