求助，伞状采样画自由能曲线，最后曲线急剧升高，不平衡，怎么处理？

sakura

直方图看过了，没什么问题。

霜晨月

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.

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

k64_cc

霜晨月 发表于 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如果是的话，按理说结果应该没啥大问题……不过还是建议并行做几组看看重现性。

霜晨月

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相比，形状倒是有点相似，但纵坐标和极大值极小值的横坐标就差远了。


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

k64_cc

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

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

文件名 center bias deltaT

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

霜晨月

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

http://membrane.urm ...

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

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

k64_cc

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

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

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

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

霜晨月

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

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

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

k64_cc

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

跑不少于3组，算误差。

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

霜晨月

借楼请教下，用伞形采样算PMF，做error bar是必需的吗？error bar该怎么做，是不是要平行跑2~3个repeat（设定不同的初始速度种子数）？@k64_cc @fhh2626


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

wasngsimin

我想问一下计算用伞状抽样+WHAM计算PMF时需要分子动力学软件提供给WHAM什么信息呢？

sakura

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

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

sakura

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

曲线不好，就没算偏差。

k64_cc

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

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

k64_cc

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

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

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