求助：伞形采样结果分析，对距离等数据求构象平均值有意义吗？

霜晨月

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

k64_cc

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

霜晨月

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

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

fhh2626

1、对整个轨迹中的采样做加权
2、对每个bin中的采样做平均

这两个都是常规分析

霜晨月

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

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

fhh2626

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

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

k64_cc

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

MBAR，MBAR……

霜晨月

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

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

霜晨月

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

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

k64_cc

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

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

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

霜晨月

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中的相应函数代码，发现完全看不懂，太难了。。。

k64_cc

霜晨月 发表于 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。

霜晨月

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

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

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

k64_cc

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

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

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

霜晨月

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

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