计算化学公社

标题: 使用MULE在势能面/自由能面上寻找最低反应路径 [打印本页]

作者
Author: fhh2626 时间: 2020-6-3 22:28
标题: 使用MULE在势能面/自由能面上寻找最低反应路径
本帖最后由 fhh2626 于 2020-6-3 22:34 编辑

MUltidimensional Least Energy pathway (MULE)的原理其实就是Dijkstra算法，属于贪心算法的一种，具体数学内容就不多扯了。和其他软件相比，MULE有几个优点：
• 得到的一定是理论严谨的最低自由能路径，而其他方法，比如LFEP和String会陷入局部极小值
• 可以处理周期性反应坐标（比如二面角）
• 可以处理任意维度，现在网上能找到的大多数代码都只能处理二维路径
• 可以施加约束，从而实现A*算法，或者寻找第2/3/N条最低自由能路径，亦或者寻找经过某个点的最低自由能路径

MULE的windows主程序/源代码/示例文件都在压缩包里面提供，linux和mac用户需要自己编译mule.cpp文件

使用MULE需要写一个config文件，然后执行
./mule.exe config
软件就会自动生成一个.traj和一个.energy文件，分别对应路径和沿路径的能量/自由能变化

/example1文件夹中给出了config文件的写法
[mule]
directory    = nanma_ref.pmf
initial       = -156, 160
end          =    78, -58
pbc          =    1, 1
软件会读取执行目录的nanma_ref.pmf文件，因为这个文件是NAMD生成的，软件会自己读取里面记录的上下界和格点精度等信息，如果给出的文件是普通格式的话需要人为在config文件中提供这些信息。剩下的三行分别表示路径从(-156, 160)开始，到(78, -58)为止，反应坐标的两个维度都是周期性的，如果非周期性的反应坐标则将pbc设为0。找出来的路径如图中左上所示。

在/example2文件夹中给出了找这个文件中第2低的反应路径的例子，config文件内容如下：
[mule]
directory    = nanma_ref.pmf
initial       = -156, 160
end          =    78, -58
lowerboundary = -180, -180
upperboundary = 178, 178
width          =    2, 2
pbc          =    1, 1
writeExploredPoints = 1
target       =    136, -32, -0.1, -0.1
这个例子中我提供了一个普通dat格式的文件，所以需要手动设置上下界和格点精度(lowerboundary, upperboundary and width)，writeExploredPoints选项设为1会让程序输出探索过的所有点，一般来说设计A*算法时有些用，这里没用，不去管它。值得解释的是最后一行，这一行施加了一个Manhattan约束，具体数学公式不多说了。这里代表采用(0.1, 0.1)的力常数将探索过程推离(136, -32)，其中(136, -32)是example1中得到的过渡态。这个约束会使得算法找到一条较为远离上一步中过渡态的路径，如右上子图所示

如果我们看这个自由能图，第一反应的话一定是找到从中间走的这个路径，example3给出了寻找这条路径的方法
[mule]
directory    = nanma_ref.pmf
initial       = -156, 160
end          =    78, -58
lowerboundary = -180, -180
upperboundary = 178, 178
width          =    2, 2
pbc          =    1, 1
writeExploredPoints = 1
target       =    0, 100, 0.1, 0.1
其他参数都一样，最后一行代表采用(0.1, 0.1)的力常数将探索过程拉向(0, 100)这个点（接近于这条路径的过渡态），得到的路径为左下图

右下的图对比了三条路径的自由能变化，可以看见第一条路径明显最有利
(, 下载次数 Times of downloads: 185)

Reference:       Fu H., Chen H, Wang X, et al. Finding an Optimal Pathway on a Multidimensional Free-Energy Landscape. J. Chem. Inf. Model. 2020. doi:10.1021/acs.jcim.0c00279

(, 下载次数 Times of downloads: 303)

作者
Author: k64_cc 时间: 2021-2-15 20:10
我说这项目为什么看着如此眼熟，原来三作是我本科室友……

作者
Author: youyno 时间: 2023-2-10 17:45
三维势能面可以用这个程序找到极小值点并搜索到连接两个点之间的路径吗？

作者
Author: fhh2626 时间: 2023-2-13 11:38

youyno 发表于 2023-2-10 17:45
三维势能面可以用这个程序找到极小值点并搜索到连接两个点之间的路径吗？

可以

作者
Author: elpa 时间: 2023-2-15 11:31
请问如何在寻找路径前确定极小点位置？

作者
Author: youyno 时间: 2023-2-15 19:54
本帖最后由 youyno 于 2023-2-15 19:59 编辑

elpa 发表于 2023-2-15 11:31
请问如何在寻找路径前确定极小点位置？

同问，pmf文件是不是就是坐标和自由能信息？

作者
Author: youyno 时间: 2023-2-15 20:08
将plumed生成的自由能文件编辑成pmf文件，试了下还是不行，报This is not an NAMD PMF file!

作者
Author: fhh2626 时间: 2023-2-16 10:37

youyno 发表于 2023-2-15 20:08
将plumed生成的自由能文件编辑成pmf文件，试了下还是不行，报This is not an NAMD PMF file!

看example2。如果是NAMD文件程序会自动读取边界等信息，如果不是的话需要在config文件里面手动提供

作者
Author: elpa 时间: 2023-2-16 15:52

youyno 发表于 2023-2-15 20:08
将plumed生成的自由能文件编辑成pmf文件，试了下还是不行，报This is not an NAMD PMF file!

我也是plumed生成的fes.dat自由能文件，问题是如何转换成pmf文件？

作者
Author: youyno 时间: 2023-2-16 16:46

elpa 发表于 2023-2-16 15:52
我也是plumed生成的fes.dat自由能文件，问题是如何转换成pmf文件？

把最后几列不需要的删掉后再改个后缀就行，关键是怎么确定极小值点

作者
Author: elpa 时间: 2023-2-16 19:53

youyno 发表于 2023-2-16 16:46
把最后几列不需要的删掉后再改个后缀就行，关键是怎么确定极小值点

我用这个http://www.metadynamics.cz/metadynminer/index.html#par3.1 找的极小值，函数fesminima() 可以。

作者
Author: elpa 时间: 2023-2-16 19:58

youyno 发表于 2023-2-16 16:46
把最后几列不需要的删掉后再改个后缀就行，关键是怎么确定极小值点

最后几列具体是哪几列？应该一共5列，最后两列是微分。是不是只留前三列？能不能给我一个参考示例，fes.dat转pmf。

作者
Author: mz121 时间: 2024-1-22 20:50
您好，我的自由能并不是填满整个空间，mule在寻找图中箭头起点到终点的路径时，会把没有数据得空白部分认为是0，得到经过空白部分的路径，但实际上这里是没有数据的，这该如何解决呢？

作者
Author: fhh2626 时间: 2024-1-26 09:36

mz121 发表于 2024-1-22 20:50
您好，我的自由能并不是填满整个空间，mule在寻找图中箭头起点到终点的路径时，会把没有数据得空白部分认为 ...

你写个脚本，让自由能面填满整个空间，没有值的地方都填个很大的数就可以（Colvars和Plumed输出的自由能面应该都是有值的，你这个怎么弄出来的）。

作者
Author: fhh2626 时间: 2024-1-26 09:38

youyno 发表于 2023-2-16 16:46
把最后几列不需要的删掉后再改个后缀就行，关键是怎么确定极小值点

可以写个简单的贪心算法程序找，另外即使不选极小值点，MULE找出来的路径肯定也会经过附近的极小值点，直接从输出文件中找就行

作者
Author: fhh2626 时间: 2024-1-26 09:39

elpa 发表于 2023-2-16 19:58
最后几列具体是哪几列？应该一共5列，最后两列是微分。是不是只留前三列？能不能给我一个参考示例，fes.d ...

PMF和dat是一个格式，就是

CV1 CV2 ... CVn deltaG

作者
Author: mz121 时间: 2024-1-26 13:27

fhh2626 发表于 2024-1-26 09:36
你写个脚本，让自由能面填满整个空间，没有值的地方都填个很大的数就可以（Colvars和Plumed输出的自由能 ...

这个自由能是通过马尔可夫链得出的。我写了个脚本，填充了空间，但mule查找出来的自由能路径上的自由能值为0。我自由能文件中并未有0值，不知道为何mule会得到0值？

作者
Author: mz121 时间: 2024-1-26 13:35
本帖最后由 mz121 于 2024-1-26 13:36 编辑

fhh2626 发表于 2024-1-26 09:36
你写个脚本，让自由能面填满整个空间，没有值的地方都填个很大的数就可以（Colvars和Plumed输出的自由能 ...

这个自由能面是通过构建马尔科夫链得到的，现在我已经填满了空间，但mule查找的自由能路径上的自由能值大多数为0，我自由能面中没有0值，不知道mule为何会查找出0值？

作者
Author: fhh2626 时间: 2024-1-28 11:44

mz121 发表于 2024-1-26 13:35
这个自由能面是通过构建马尔科夫链得到的，现在我已经填满了空间，但mule查找的自由能路径上的自由能值大 ...

你的dat/pmf文件是等间隔的吗，你看看例子中的文件，弄成一样的格式应该就行

作者
Author: gehan 时间: 2024-4-16 10:12
请问楼主，左上，左下和右上三个图是怎么绘制出的，二维等值面可以使用origin生成，图中的黑色路径的线是怎么在图中显示出来的，还有就右下的自由能变化曲线的横坐标是取什么为变量的呢？谢谢

作者
Author: fhh2626 时间: 2024-4-16 14:22

gehan 发表于 2024-4-16 10:12
请问楼主，左上，左下和右上三个图是怎么绘制出的，二维等值面可以使用origin生成，图中的黑色路径的线是怎 ...

origin可以添加图层，添加一个图层画散点图就行了

右下角的图是以路径为坐标，如果你用的MULE的话会直接输出这个文件

作者
Author: gehan 时间: 2024-4-18 10:01

fhh2626 发表于 2024-4-16 14:22
origin可以添加图层，添加一个图层画散点图就行了

右下角的图是以路径为坐标，如果你用的MULE的话会直 ...

谢谢您的回答，但是我使用MULE生成了.energy文件里是只有一列能量的数值，没有横坐标哎，不知道怎么设置路径为坐标。还有一个问题，请问有没有什么办法，找到势能面中的极小值点以及过渡态的位置后，怎么才能找到对应的结构呢？

作者
Author: fhh2626 时间: 2024-4-18 10:15

gehan 发表于 2024-4-18 10:01
谢谢您的回答，但是我使用MULE生成了.energy文件里是只有一列能量的数值，没有横坐标哎，不知道怎么设置 ...

横坐标就是0-1自己插值

MD会输出轨迹文件啊，在轨迹文件里面找

作者
Author: gehan 时间: 2024-4-18 10:48

fhh2626 发表于 2024-4-18 10:15
横坐标就是0-1自己插值

MD会输出轨迹文件啊，在轨迹文件里面找

再次感谢，那请问是怎么根据势能面上极值点的信息，定位到MD轨迹文件中的哪一帧结构呢，是根据极值点的对应的CV1和CV2，然后去fes.dat,或者该文章中的.pmf文件中找到对应的值，然后在该文件中找到对应的帧数吗。还是说是自己看MD轨迹动画去找到想要极值点结构呢。

作者
Author: fhh2626 时间: 2024-4-18 11:54

gehan 发表于 2024-4-18 10:48
再次感谢，那请问是怎么根据势能面上极值点的信息，定位到MD轨迹文件中的哪一帧结构呢，是根据极值点的对 ...

你知道极值点就知道CVs了，就可以去轨迹里筛选结构

作者
Author: zhonghaoxuan 时间: 2025-6-23 10:48
老师您好，跟您请教一下，使用MULE的时候，我生成的反应路径上只有70个数据点，能量也只有70个，请问如何设置CONFIG文件可以使数据点更多呢

作者
Author: fhh2626 时间: 2025-6-27 20:06

zhonghaoxuan 发表于 2025-6-23 10:48
老师您好，跟您请教一下，使用MULE的时候，我生成的反应路径上只有70个数据点，能量也只有70个，请问如何设 ...

这是你提供的自由能面的精度决定的，你可以做一个插值或者平滑

作者
Author: zhonghaoxuan 时间: 2025-7-3 09:37

fhh2626 发表于 2025-6-27 20:06
这是你提供的自由能面的精度决定的，你可以做一个插值或者平滑

非常感谢

欢迎光临计算化学公社 (http://bbs.keinsci.com/)