计算化学公社

标题: 求助：伞形采样结果分析，对距离等数据求构象平均值有意义吗？ [打印本页]

作者
Author: 霜晨月 时间: 2022-9-9 12:59
标题: 求助：伞形采样结果分析，对距离等数据求构象平均值有意义吗？
各位老师，我做了个蛋白质的伞形采样，审稿人要求分析distance、dihedral等数据的构象平均值以及fluctuations。但伞形采样是biased simulation，每个帧的偏置力、权重等并不一样，这样对每帧求构象平均值和fluctuations有意义吗？
此外，审稿人还要求作出各个数据（distance、dihedral等）对整个反应坐标的变化曲线，但伞形采样是分成很多窗口来进行模拟，每个窗口只cover反应坐标的一小段，并没有完整的从Reactant到Product的模拟过程。这个该怎么做？能否把所有窗口的轨迹合并，将每一帧的distances等参数对整个反应坐标做散点图？（但这样似乎又遇到上面的问题，就是每一帧的权重都不一样。）

作者
Author: k64_cc 时间: 2022-9-9 13:51
本帖最后由 k64_cc 于 2022-9-9 13:53 编辑

有个神奇的东西叫做MBAR，可以把不同bias下的所有轨迹纳入同一个统计，解决上述所有问题。
但是如果他想看到在反应坐标取特定值的时候其他结构信息的分布，那最好的选择其实是加非常强的restraint重跑一遍。

作者
Author: 霜晨月 时间: 2022-9-9 14:08

k64_cc 发表于 2022-9-9 13:51
有个神奇的东西叫做MBAR，可以把不同bias下的所有轨迹纳入同一个统计，解决上述所有问题。
但是如果他想看 ...

我理解审稿人的意思大概是，canonical MD文章都有这些分析，所以umbrella sampling也要有。。。
倒是见过一些文章，对umbrella sampling也做了这些图，但似乎没提数据统计处理的问题，也许是直接把各帧的数据拿来作图了。所以不知道这样作图有无意义

作者
Author: fhh2626 时间: 2022-9-9 14:27
1、对整个轨迹中的采样做加权
2、对每个bin中的采样做平均

这两个都是常规分析

作者
Author: 霜晨月 时间: 2022-9-9 14:47

fhh2626 发表于 2022-9-9 14:27
1、对整个轨迹中的采样做加权
2、对每个bin中的采样做平均

每个bin中的采样可以直接平均？不需要加权吗？
其实我就是不知道如何对每个构象进行加权

作者
Author: fhh2626 时间: 2022-9-9 14:50

霜晨月发表于 2022-9-9 14:47
每个bin中的采样可以直接平均？不需要加权吗？
其实我就是不知道如何对每个构象进行加权

每个bin对应的自由能值（或者施加的偏置力）应该是一样的，也就是说权重是相同的，当然如果是两个窗口的重叠部分的话则需要各自计算

作者
Author: k64_cc 时间: 2022-9-9 14:52

霜晨月发表于 2022-9-9 14:47
每个bin中的采样可以直接平均？不需要加权吗？
其实我就是不知道如何对每个构象进行加权

MBAR，MBAR……

作者
Author: 霜晨月 时间: 2022-9-9 15:18

k64_cc 发表于 2022-9-9 14:52
MBAR，MBAR……

我之前研究过MBAR，这是个python库，自带的脚本只有计算PMF的，要想reweight需要自己写脚本
看来得专门学习一段时间的python了

作者
Author: 霜晨月 时间: 2022-9-9 15:42

fhh2626 发表于 2022-9-9 14:50
每个bin对应的自由能值（或者施加的偏置力）应该是一样的，也就是说权重是相同的，当然如果是两个窗口的 ...

比如同一个bin里面有两个构象，构象A是在窗口1的中央，构象B是在窗口1的边缘，A和B的偏置力和权重肯定是不一样的吧？

作者
Author: k64_cc 时间: 2022-9-9 19:07

霜晨月发表于 2022-9-9 15:18
我之前研究过MBAR，这是个python库，自带的脚本只有计算PMF的，要想reweight需要自己写脚本
看来得专门 ...

MBAR是个算法，你看的那个叫pymbar。

用MBAR处理US的PMF，本质上就是把每个state放一起做reweighting，然后给出一个没有bias的估计。他给出无偏的PMF的时候，已经算出来你需要的weight了。

作者
Author: 霜晨月 时间: 2022-9-9 22:49

k64_cc 发表于 2022-9-9 19:07
MBAR是个算法，你看的那个叫pymbar。

用MBAR处理US的PMF，本质上就是把每个state放一起做reweighting ...

pymbar调用MBAR计算PMF的输出结果如下：

Binning data...
Evaluating reduced potential energies...
Running MBAR...
K (total states) = 17, total samples = 142714
N_k =
[11314  9990 11112  8861  7883  4661  7781  8741  8105  8059  7063  8092
  8317  7566  8936  8067  8166]
There are 17 states with samples.
Initializing free energies to zero.
Initial dimensionless free energies with method zeros
f_k =
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
Determining dimensionless free energies by Newton-Raphson / self-consistent iteration.
self consistent iteration gradient norm is 9.0031e+07, Newton-Raphson gradient norm is 6.3272e+08
Choosing self-consistent iteration on iteration 0
self consistent iteration gradient norm is 5.7465e+07, Newton-Raphson gradient norm is 1.8677e+07
Choosing self-consistent iteration for lower gradient on iteration 1
self consistent iteration gradient norm is 3.9352e+07, Newton-Raphson gradient norm is 2.2391e+06
Newton-Raphson used on iteration 2
self consistent iteration gradient norm is 3.8027e+05, Newton-Raphson gradient norm is    787.14
Newton-Raphson used on iteration 3
self consistent iteration gradient norm is    101.9, Newton-Raphson gradient norm is 4.4636e-05
Newton-Raphson used on iteration 4
self consistent iteration gradient norm is 6.5189e-06, Newton-Raphson gradient norm is 1.6842e-19
Newton-Raphson used on iteration 5
self consistent iteration gradient norm is 2.3465e-20, Newton-Raphson gradient norm is 1.9608e-22
Newton-Raphson used on iteration 6
Converged to tolerance of 1.257696e-14 in 7 iterations.
Of 7 iterations, 5 were Newton-Raphson iterations and 2 were self-consistent iterations
Final dimensionless free energies
f_k =
[  0.       -6.35575236 -10.35715215 -12.5221046  -13.38232867
-13.51146418 -13.41177059 -13.78373329 -14.5068132  -15.38831936
-16.27500615 -17.20686221 -18.04267401 -18.68934648 -19.12060088
-19.19076083 -19.03618277]
MBAR initialization complete.
PMF (in units of kT)
   bin       f    df
1.81 24.835 0.512
1.83 24.562 0.461
1.84 23.407 0.290
…………（下略）…………
请问老师，哪部分是reweight的结果呢？是红字部分吗？但这个结果可以解决主楼的问题吗？
（只有红色部分看上去像是reweight结果，但刚刚检查了一下，其实这是17个窗口的帧数除以相应窗口的correlation time得到的“有效帧数”，但有几个窗口的数值对不上）

刚才追查了mbar.py中的相应函数代码，发现完全看不懂，太难了。。。

作者
Author: k64_cc 时间: 2022-9-10 07:33

霜晨月发表于 2022-9-9 22:49
pymbar调用MBAR计算PMF的输出结果如下：

Binning data...

你做u_kln的时候多加一个state，是不带bias的无偏state，有效样本数是0，用这个u_kln矩阵创建MBAR实例。MBAR迭代收敛后，其中的W_nk attribute就是你要的weight。

或者你也可以用你贴出来的script算PMF，然后用MBAR.compute_multiple_expectations这个函数，把每一个frame的histogram做reweight。

作者
Author: 霜晨月 时间: 2022-9-10 11:44

k64_cc 发表于 2022-9-10 07:33
你做u_kln的时候多加一个state，是不带bias的无偏state，有效样本数是0，用这个u_kln矩阵创建MBAR实例。M ...

谢谢老师，我再去研究下代码

有效样本数怎么可能是0？除非用一个空文件？

作者
Author: k64_cc 时间: 2022-9-10 15:24

霜晨月发表于 2022-9-10 11:44
谢谢老师，我再去研究下代码

有效样本数怎么可能是0？除非用一个空文件？

不对它采样就是0了。MBAR输入文件的能量矩阵其实是要你自己构造的，你用的script实现了数据处理的部分，但是对于这种没在他们例子里的需求，数据处理显然要自己重新实现。

作者
Author: 霜晨月 时间: 2022-9-10 16:19

k64_cc 发表于 2022-9-10 15:24
不对它采样就是0了。MBAR输入文件的能量矩阵其实是要你自己构造的，你用的script实现了数据处理的部分， ...

谢谢老师，我再仔细研究下

作者
Author: fhh2626 时间: 2022-9-10 17:15

霜晨月发表于 2022-9-9 15:42
比如同一个bin里面有两个构象，构象A是在窗口1的中央，构象B是在窗口1的边缘，A和B的偏置力和权重肯定是 ...

bin不是窗口，bin是窗口的再细分，你可以把bin分得无限细，但是你画出来的图的精度总是有限的

你如果每一帧的偏置力或者反应坐标值都存下来的话你也可以手动加权

作者
Author: 霜晨月 时间: 2022-9-17 13:20
本帖最后由霜晨月于 2022-9-17 13:21 编辑

fhh2626 发表于 2022-9-10 17:15
bin不是窗口，bin是窗口的再细分，你可以把bin分得无限细，但是你画出来的图的精度总是有限的

你如果 ...

老师好，我把MBAR的W_nk提取出来，得到了每一帧的weight。接下来我想这样做：

1. 将整个反应坐标均匀分割成很多bin
2. 将各个帧（各自带有权重）按照CV值assign到各个bin中
3. 将每个bin里面的所有帧的距离等数值按照权重取平均值，并计算fluctuation
4. 将整个反应坐标上每个bin的平均距离值和fluctuation描绘曲线

请问第2步，将各个帧按照其CV值assign到各个bin中，这样做对不对？因为我记得用wham等方法计算PMF时，各个帧不一定是按照其CV值归属到bin的，比如某帧CV=0.85，不一定属于0.84~0.86的bin，而可能被assign到其他的bin。（这一点我记忆不确，是看的一篇中文介绍资料，现在看具体的PMF计算代码看不懂了。）

@k64_cc

作者
Author: fhh2626 时间: 2022-9-20 20:26

霜晨月发表于 2022-9-17 13:20
老师好，我把MBAR的W_nk提取出来，得到了每一帧的weight。接下来我想这样做：

1. 将整个反应坐标均匀 ...

bin就是按照CV来分的，不过如果不同窗口之间有overlap的话，加权的时候要妥善处理重叠部分

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