葫芦脲-质子化己二胺的解离自由能，伞状采样计算PMF后，求直方图不连续的解决方法

JCenter

请教论坛的各位大佬，我最近在计算水中葫芦脲-质子化己二胺的解离自由能。流程如下：我将体系在298k和1bar下进行1ns预平衡，然后通过拉伸分子动力学，将质子化己二胺从葫芦脲空腔中牵引出来，两者质心距离增加至2nm，如图1（水分子没有显示），牵引过程的参数文件为MD_pull.mdp，伞形牵引的简谐力常数为600 kJ/(mol*nm2)。牵引完后，每隔约0.1nm作为一次采样，得到了不同质心距离下的初始构型。我将每个窗口下的初始构型均进行500psNPT预平衡和10ns的MD模拟。NPT预平衡的动力学参数文件为NPT_umbrella.mdp，MD模拟的动力学参数文件为MD_umbrella.mdp,每个窗口下pull_coord1_rate均改为0，pull_coord1_k还是600 kJ/(mol*nm2)。然而，当我计算完成后，使用gmx wham -it tpr_files.dat -if pullf_files.dat -b 500 -o -hist（舍弃MD模拟的前500ps）和qtgrace -nxy histo.xvg，发现直方图不连续（0.15-0.5nm），而且运行gmx wham时产生了一些warning,就是直方图中缺少的部分，图2。为了其连续性，对于两组质心距离在0.15-0.5nm的几个初始构型，我将参数文件中的pull_coord1_k从600均改为10000 kJ/(mol*nm2)，重新跑了NPT预平衡和MD,结果直方图依旧不连续，如图3。不连续的直方图得到的PMF曲线在0.2-0.5nm不连续，图4。这个现象如何有效解决，使得直方图中各窗口均相互重叠和PMF曲线连续呢，还是说，这样的数据也可以使用呢。感谢大家！（感觉直方图中缺失的部分应该是质子化己二胺的质子化胺基与葫芦脲端口的6个羰基相互作用的结果）

fhh2626

首先，你这个自由能计算是收敛不了的，你选择的CV(s)无法充分描述你要研究的运动。要记住自由能计算是可逆的，你选的反应坐标要能驱动你研究的运动过程的可逆发生。这个体系如果你只考虑质心距，不考虑其他球面坐标的话，如果你用MtD或者ABF，就会发现你只能发现配体被拉出来，拉不回去了。

其次，这个过程很多年前就被研究过了
https://pubs.rsc.org/en/content/articlelanding/2014/cp/c4cp04200j

JCenter

fhh2626 发表于 2024-4-26 09:51
首先，你这个自由能计算是收敛不了的，你选择的CV(s)无法充分描述你要研究的运动。要记住自由能计算是可逆 ...

好的，感谢大佬，文献我好好学习一下。此外，您讲的内容里，学生有几点疑惑，希望大佬您有机会指点一下。1.“你这个自由能计算是收敛不了的，你选择的CV(s)无法充分描述你要研究的运动。”这个自由能计算的收敛性是如何判断的呢。2。关于点群变量CV(s),和考虑的反应坐标。如何考虑球面坐标以充分描述我体系的运动，是选择pull代码中pull_coord1_geometry哪种方式呢，通过看代码，发现没有合适的关键词。

fhh2626

JCenter 发表于 2024-4-26 11:39
好的，感谢大佬，文献我好好学习一下。此外，您讲的内容里，学生有几点疑惑，希望大佬您有机会指点 ...

这个一两句话说不清楚啊，你可以先学一下专业一点的自由能计算模块，比如Colvars和Plumed，都可以和GMX联用，里面有eABF或者MtD，并且可以更加复杂的甚至自定义的CVs。如果做得糙一些，也可以施加几何约束

JCenter

fhh2626 发表于 2024-4-26 11:51
这个一两句话说不清楚啊，你可以先学一下专业一点的自由能计算模块，比如Colvars和Plumed，都可以和GMX联 ...

好嘞，感谢大佬

七尺贱

fhh2626 发表于 2024-4-26 11:51
这个一两句话说不清楚啊，你可以先学一下专业一点的自由能计算模块，比如Colvars和Plumed，都可以和GMX联 ...

Plumed跟着手册学吗，有没有其它学习教程推荐

fhh2626

七尺贱 发表于 2024-4-26 18:29
Plumed跟着手册学吗，有没有其它学习教程推荐

自带的教程挺详细了吧

JCenter

fhh2626 发表于 2024-4-26 09:51
首先，你这个自由能计算是收敛不了的，你选择的CV(s)无法充分描述你要研究的运动。要记住自由能计算是可逆 ...

大佬，想请教下您：我最近重复了一下US方法，先简单的使用质心距离作为反应坐标，不过是将葫芦脲进行了位置限制。反应坐标每隔0.1nm进行一次采样。每个采样窗口均将环状分子葫芦脲进行了位置限制。弹簧力常数为1500，每个窗口均进行了1ns的npt预平衡和20 ns的md产生相。然而通过wham得到的直方图仍然显示缺失，如图所示。是US方法本身的问题吗，还是说我的某些参数没有使得采样充分呢。我如果后面用插件，使用CVs包括葫芦脲的球面坐标，感觉还是会碰到采样缺失的问题。如果说是US方法本身的问题，是不是eABF会好一点呢

pixy

你好，我也遇到了类似的问题，想请教你这个问题是怎么解决的？非常感谢！

JCenter

pixy 发表于 2024-9-6 03:08
你好，我也遇到了类似的问题，想请教你这个问题是怎么解决的？非常感谢！

高能垒的地方由于力常数过小，导致该处没法被拴住，采样不充分。因此要在高能垒处的窗口采用较高的K值。

pixy

JCenter 发表于 2024-9-9 15:31
高能垒的地方由于力常数过小，导致该处没法被拴住，采样不充分。因此要在高能垒处的窗口采用较高的K值。

嗯嗯是的，我最近也在尝试你说的这种方法，但是我不确定力常数太大的话会不会影响能垒的准确性？这个K值一般有上限吗？设5000左右会不会太高呢？

JCenter

pixy 发表于 2024-9-9 22:44
嗯嗯是的，我最近也在尝试你说的这种方法，但是我不确定力常数太大的话会不会影响能垒的准确性？这个K值 ...

不会的，我都有设置过k=15000的情况

pixy

JCenter 发表于 2024-9-10 00:40
不会的，我都有设置过k=15000的情况

嗯嗯，好的，非常感谢~

ljh123

JCenter 发表于 2024-9-9 15:31
高能垒的地方由于力常数过小，导致该处没法被拴住，采样不充分。因此要在高能垒处的窗口采用较高的K值。

请问不同窗口可以用不同的力常数吗