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

霜晨月

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

http://membrane.urm ...

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

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

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-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 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-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.

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