计算化学公社

标题: 求助，伞状采样画自由能曲线，最后曲线急剧升高，不平衡，怎么处理？ [打印本页]

作者
Author: sakura 时间: 2019-2-20 15:16
标题: 求助，伞状采样画自由能曲线，最后曲线急剧升高，不平衡，怎么处理？
直方图看过了，没什么问题。

作者
Author: fhh2626 时间: 2019-2-20 17:28
自由能曲线就是直方图算出来的，直方图和自由能曲线是一回事，直方图“没问题”自由能曲线就该“没问题”

作者
Author: sakura 时间: 2019-2-20 19:38

fhh2626 发表于 2019-2-20 17:28
自由能曲线就是直方图算出来的，直方图和自由能曲线是一回事，直方图“没问题”自由能曲线就该“没问题”

那这条曲线怎么解释？

作者
Author: k64_cc 时间: 2019-2-21 11:34
首先自己说没问题不一定就是真没问题。最好把overlapping发上来看一眼，万一你不知道什么是有问题，岂不是白费事。

你这个问题的第一种可能是后面数据不足，这个没什么办法，加window加bias吧。
第二种可能是没well sampled，你把thermostat换langevin，再加数据，看能不能降下来。
第三种可能是它就应该是这样，建议回家反思人品。

作者
Author: sakura 时间: 2019-2-21 14:32

k64_cc 发表于 2019-2-21 11:34
首先自己说没问题不一定就是真没问题。最好把overlapping发上来看一眼，万一你不知道什么是有问题，岂不是 ...

谢谢，我计算该晶体另一个面的自由能的时候，曲线平衡的很好。
到了这个面，相同的条件，就变成这个样子了。

作者
Author: hdmd 时间: 2019-2-21 15:29
通过这两个图可以看出，转化过程没问题（至少在你给出的部分），是胶束分离还是跨越孔道？或者其他过程，反应坐标是什么，采样值是什么，偏置力如何设置的，有没有可能这个过程就是这样的？

作者
Author: k64_cc 时间: 2019-2-21 15:36

sakura 发表于 2019-2-21 14:32
谢谢，我计算该晶体另一个面的自由能的时候，曲线平衡的很好。
到了这个面，相同的条件，就变成这个样子 ...

你有没有觉得2.0左右贼不均匀…

要是峰的shift比较大，建议加bias让它别太偏。虽然理论上没影响，但是峰偏离bias中心点太多的话，需要的数据量是真的多…

作者
Author: k64_cc 时间: 2019-2-21 15:37

sakura 发表于 2019-2-21 14:32
谢谢，我计算该晶体另一个面的自由能的时候，曲线平衡的很好。
到了这个面，相同的条件，就变成这个样子 ...

话说你PMF倒是做error bar啊……

作者
Author: sakura 时间: 2019-2-21 16:09

k64_cc 发表于 2019-2-21 15:37
话说你PMF倒是做error bar啊……

曲线不好，就没算偏差。

作者
Author: sakura 时间: 2019-2-21 16:11

hdmd 发表于 2019-2-21 15:29
通过这两个图可以看出，转化过程没问题（至少在你给出的部分），是胶束分离还是跨越孔道？或者其他过程， ...

是，从晶体表面拉出一个离子，计算自由能。这是这个晶体另一个表面的自由能曲线

作者
Author: wasngsimin 时间: 2019-2-21 21:06
我想问一下计算用伞状抽样+WHAM计算PMF时需要分子动力学软件提供给WHAM什么信息呢？

作者
Author: 霜晨月 时间: 2020-1-13 20:08
借楼请教下，用伞形采样算PMF，做error bar是必需的吗？error bar该怎么做，是不是要平行跑2~3个repeat（设定不同的初始速度种子数）？@k64_cc @fhh2626

另外，Amber里面怎么做5楼的这种overlapping图，这个应该是各个窗口的概率密度曲线重叠在一起，gmx wham可以自动给出这个，Amber里面没找到这种工具，是不是得自己写个小程序

作者
Author: k64_cc 时间: 2020-1-15 11:53

霜晨月发表于 2020-1-13 20:08
借楼请教下，用伞形采样算PMF，做error bar是必需的吗？error bar该怎么做，是不是要平行跑2~3个repeat（设 ...

跑不少于3组，算误差。

做overlapping确实得自己写……读CV，画直方图，把所有window的数据都画一起。

作者
Author: 霜晨月 时间: 2020-1-15 15:23

k64_cc 发表于 2020-1-15 11:53
跑不少于3组，算误差。

做overlapping确实得自己写……读CV，画直方图，把所有window的数据都画一起。

哦哦，谢谢。
这样说来工作量翻了三倍啊

作者
Author: k64_cc 时间: 2020-1-15 16:55

霜晨月发表于 2020-1-15 15:23
哦哦，谢谢。
这样说来工作量翻了三倍啊

2020年了，弄个显卡，好像还行？或者也可以bootstrap，但是需要的总数据量不变。

http://membrane.urmc.rochester.edu/?page_id=126

他们的wham工具带了这个功能，具体看文档

作者
Author: 霜晨月 时间: 2020-1-16 03:19

k64_cc 发表于 2020-1-15 16:55
2020年了，弄个显卡，好像还行？或者也可以bootstrap，但是需要的总数据量不变。

http://membrane.urm ...

是QM/MM MD，显卡用不上。。。（如果不是要用QM，我其实更习惯用gmx）
bootstrap不熟悉，打算研究一下。

我现在用的就是他的这个wham，输出有+/-这项，但都是零，是因为没做repeat吧？

作者
Author: k64_cc 时间: 2020-1-16 11:04
本帖最后由 k64_cc 于 2020-1-16 11:05 编辑

霜晨月发表于 2020-1-16 03:19
是QM/MM MD，显卡用不上。。。（如果不是要用QM，我其实更习惯用gmx）
bootstrap不熟悉，打算研究一下。 ...

因为你whamlist少了一项。格式应该是：

文件名 center bias deltaT

deltaT那项是说你的体系的轨迹经过多少个frame之后与第一个frame不相关。你加上那项，他就可以把轨迹按这个分段，然后重新抽样。
不过QMMM的话……您还是用三个不同初始态跑仨轨迹吧。

作者
Author: 霜晨月 时间: 2020-1-16 17:50
本帖最后由霜晨月于 2020-1-16 17:52 编辑

k64_cc 发表于 2020-1-16 11:04
因为你whamlist少了一项。格式应该是：

文件名 center bias deltaT

这个“你的体系的轨迹经过多少个frame之后与第一个frame不相关”，是不是有些文献上的correlation time啊？

我这两天倒弄了一个工具叫Pymbar，它是先读取CV，利用一个函数timeseries.statisticalInefficiency(CV)来算出correlation time，然后给出PMF曲线和error bar。不过它算出来的correlation time都很小，一般都是一点几帧，最大也就是五点几帧。（感觉好像不是我的数据问题，因为它自带的example算出来结果就是如此。）

这个是我用Pymbar做出来的PMF和error bar，20个窗口，每个窗口120ps，80个bin。与Grossfield的WHAM给出的PMF相比，形状倒是有点相似，但纵坐标和极大值极小值的横坐标就差远了。
(, 下载次数 Times of downloads: 30)

上面是Pymbar，下面是WHAM，纵坐标都是kcal/mol（Pymbar给出的单位是kBT，已经换算了）。是不是Pymbar的结果不合理？

作者
Author: k64_cc 时间: 2020-1-17 14:45
本帖最后由 k64_cc 于 2020-1-17 14:50 编辑

霜晨月发表于 2020-1-16 17:50
这个“你的体系的轨迹经过多少个frame之后与第一个frame不相关”，是不是有些文献上的correlation time啊 ...

对就是correlation time。

你是照着这个做的吗？
https://github.com/choderalab/pymbar-examples/blob/master/umbrella-sampling-pmf/umbrella-sampling.py如果是的话，按理说结果应该没啥大问题……不过还是建议并行做几组看看重现性。

作者
Author: 霜晨月 时间: 2020-1-19 14:17

k64_cc 发表于 2020-1-17 14:45
对就是correlation time。

你是照着这个做的吗？

对，就是按这个做的

例子的CV是chi torsion，我的CV是键长差，所以自己修改了代码。
改代码时，我比较没把握的地方就是correlation time的计算，因为搞不清楚这个东西的计算原理。例子是用torsion的cos和sin来求算的，我就直接用键长差代入函数求算了。

例子的代码：
# Compute correlation times for potential energy and chi
# timeseries.  If the temperatures differ, use energies to determine samples; otherwise, use the cosine of chi

if (DifferentTemperatures):
   【略】
else:
      chi_radians = chi_kn[k,0:N_k[k]]/(180.0/numpy.pi)
      g_cos = timeseries.statisticalInefficiency(numpy.cos(chi_radians))
      g_sin = timeseries.statisticalInefficiency(numpy.sin(chi_radians))
      print "g_cos = %.1f | g_sin = %.1f" % (g_cos, g_sin)
      g_k[k] = max(g_cos, g_sin)
      print "Correlation time for set %5d is %10.3f" % (k,g_k[k])
      indices = timeseries.subsampleCorrelatedData(chi_radians, g=g_k[k])
# Subsample data.

我修改之后：
if (DifferentTemperatures):
      【略】
else:
      chi_radians = chi_kn[k,0:N_k[k]] # Read the CV values of frame n of simulation k to chi_radians
      g_k[k] = timeseries.statisticalInefficiency(chi_radians) #### 这句代码不知道是否合理
      print("Correlation time for set %5d is %10.3f" % (k,g_k[k]))
      indices = timeseries.subsampleCorrelatedData(chi_radians, g=g_k[k])
# Subsample data.

请您看看这样改合理吗？多谢多谢

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