GROMACS计算MSD的原理是什么？

jimulation

http://bbs.keinsci.com/thread-6830-1-1.html Sob老师在此贴2楼回复了MSD的计算公式。

在GROMACS的Manual中，http://manual.gromacs.org/archive/5.0.5/programs/gmx-msd.html,有这样的描述:"gmx msd computes the mean square displacement (MSD) of atoms from a set of initial positions. This provides an easy way to compute the diffusion constant using the Einstein relation. The time between the reference points for the MSD calculation is set with -trestart. The diffusion constant is calculated by least squares fitting a straight line (D*t + c) through the MSD(t) from -beginfit to -endfit (note that t is time from the reference positions, not simulation time). An error estimate given, which is the difference of the diffusion coefficients obtained from fits over the two halves of the fit interval......The diffusion coefficient is determined by linear regression of the MSD, where, unlike for the normal output of D, the times are weighted according to the number of reference points, i.e. short times have a higher weight".



有个疑问想请教大家，上述标红的reference points指的是什么，和initial positions是一回事么？"short times have a higher weight"这句该怎么理解呢？
如果我计算某100个原子的MSD，轨迹一共10ns，每10 ps存一帧，那么GROMACS计算MSD的过程，是不是以t=0时刻各原子的位置为初始位置，计算100个原子在时刻的|r(t)-r(0)|2然后取平均值？




请大家指点

sobereva

reference point应当是指下图每一个箭头的起点。-trestart设的应该是每个箭头的起点之间的间隔

jimulation

sobereva 发表于 2019-7-19 22:35
reference point应当是指下图每一个箭头的起点。-trestart设的应该是每个箭头的起点之间的间隔

谢谢Sob老师回复，清晰明了！ 还有两个问题想请教您，期待回复。
1. 我的xtc文件间隔是1 ps，而trestart是默认值（10 ps），发现得到的MSD数据间隔是1 ps。这种情况下，是否浪费了许多拟合区段呢（比如MSD曲线上，t=1 ps时刻的值，是0-1 ps，10-11 ps, 20-21 ps...取平均，而1-2 ps、2-3 ps...都没用上）？
2. 计算msd时，是不是时间t点的数量与xtc文件有关，且不可选（因为无-dt选项）。
谢谢哈！

sobereva

jimulation 发表于 2019-7-19 23:57
谢谢Sob老师回复，清晰明了！ 还有两个问题想请教您，期待回复。
1. 我的xtc文件间隔是1 ps， ...

1 可能是这样。我暂时没时间专门去看，你可以试试改小了是什么情况
2 没有-dt的话那就是如此。不设置这个选项估计是为了避免情况搞复杂

jimulation

sobereva 发表于 2019-7-20 00:06
1 可能是这样。我暂时没时间专门去看，你可以试试改小了是什么情况
2 没有-dt的话那就是如此。不设置这 ...

感谢Sob老师

cavalier

jimulation 发表于 2019-7-21 09:07
感谢Sob老师

      根据sob老师给的第二幅图，我还有一个问题始终不明白，gromacs为什么能在显示有多次restart的情况下还能画出和模拟时间一样长的msd曲线？
      我是这么理解的，要想更加精确，所选的各区段就要更短，正如图上箭头所示，就算区段最长，也只能取长度为7，亦即0-7或1-8。这样我们能画的曲线就只有t=0-7的msd曲线了，而且这种取法因为区段数量少还比较不可靠，所以应该还需再进一步增加区段数，减短每一区段长度吧。
      那么能画出来的msd曲线长度应该远远短于模拟时间才对啊？谢谢了！

sobereva

cavalier 发表于 2019-8-26 14:35
根据sob老师给的第二幅图，我还有一个问题始终不明白，gromacs为什么能在显示有多次restart的情况 ...

能给出的MSD跨度范围<= trjconv统计的轨迹长度
越大跨度的MSD采样数越少，因此越靠后越糙

cavalier

sobereva 发表于 2019-8-27 00:49
能给出的MSD跨度范围

      感谢sob老师回复！
      该不会意思是说对于一次模拟的msd曲线，需要不同msd跨度范围采样的，而后再把所有跨度的所有采样都要拿来平均？其中每种msd跨度有[（轨迹总长-跨度范围）/间隔]组采样。
    所以如果以上边第二个图：9帧，总长为8的模拟为例，是不是这样操作的：
      msd跨度范围为1的有0-1，1-2,······，7-8，这样8组采样
     msd跨度范围为2 的有0-2,1-3，······6-8，这样7组采样
      ······
      直到msd跨度为8的只有0-8，一组采样
      总计一共36组采样
      那么进行系宗平均时，36组全部都可以用来平均，来得到最终的0-1的msd曲线。而对于绘制1-2的msd曲线，能拿来平均的就只有36-8的28组了，精度也下降了。
      ······
      直到绘制7-8的msd曲线，只剩下1组，0-8的采样可以用，精度最低。
      所以最终画出来的msd曲线越往后精度越低，而刚开始的曲线又因为处在弛豫阶段也不可靠。
      好像懂了，是我说的这样吗？
      看来gromacs里那句“short times have a higher weight”我可能一直理解错了，这说的不是像0-1,1-2这种短的采样在平均的时候要有更高权重，而是在平均最终的msd曲线时越靠前的曲线，平均时使用过的采样越多。
      请问我这样理解基本正确了吗？谢谢！

sobereva

cavalier 发表于 2019-8-27 12:10
感谢sob老师回复！
      该不会意思是说对于一次模拟的msd曲线，需要不同msd跨度范围采样的，而 ...

比如计算MSD的x=1跨度的位置，相当于把0-1、1-2、2-3等共8个值求平均。x=2跨度的位置相当于把0-2、1-3、2-4等共7个值求平均。因此跨度越小，样本越多，统计误差越小。

cavalier

sobereva 发表于 2019-8-28 00:09
比如计算MSD的x=1跨度的位置，相当于把0-1、1-2、2-3等共8个值求平均。x=2跨度的位置相当于把0-2、1-3、2 ...

      感谢sob老师再次回复！
      我现在感觉“因此跨度越小，样本越多，统计误差越小。”我是可以理解的，我是认同的。
      不过，希望我没理解错。
      但是，在想要较小统计误差的前提下，不就不得不缩小跨度了。所以绘制的msd曲线肯定会小于模拟长度了？
      我想总该需要对不同跨度，进行某种共同处理、某种平均，然后整合在一起吧。比如我前面的想法，x=1至x=8的跨度都包含了0-1的msd部分，所以都可以拿来平均而得到最终的msd曲线的0-1部分。总觉得只有这样才可能做的gmx msd里模拟长度=最终输出的msd曲线长度。

sobereva

cavalier 发表于 2019-8-28 11:10
感谢sob老师再次回复！
      我现在感觉“因此跨度越小，样本越多，统计误差越小。”我是可以理 ...

“x=1至x=8的跨度都包含了0-1的msd部分”你没理解对
诸如图上x=7（即MSD=7）在计算时，利用了0-7、1-8、2-9、3-10，这里比如3-10指的是基于第3个点的结构和第10个点的结构来计算出的跨度为7的MSD值。所以并不牵扯到“0-1”

cavalier

sobereva 发表于 2019-8-28 23:19
“x=1至x=8的跨度都包含了0-1的msd部分”你没理解对
诸如图上x=7（即MSD=7）在计算时，利用了0-7、1-8、 ...

      万分感谢，我现在终于明白了，gmx msd之所以能模拟多久就输出多久的msd-t，且越往后精度越低，是因为：比如0-10的长度10的模拟，msd-t图线中，t=1是10个点平均而来，t=2是9个点平均而来······t=10则只有1个点（这个点是t=0的坐标和t=10的坐标计算出来的）。
      我发现，其实这样一来，模拟的任意两个时刻之间，其坐标都已计算过msd了，对数据的利用率已经达到百分之百了。

cavalier

sobereva 发表于 2019-8-28 23:19
“x=1至x=8的跨度都包含了0-1的msd部分”你没理解对
诸如图上x=7（即MSD=7）在计算时，利用了0-7、1-8、 ...

      sob老师，我在使用gmx msd之前其实是先自己编了一个算法，现在想来，其实我自己的算法和老师您说的差不多，但是区别，其实也是我的算法的不足之处就是数据的利用率并没有达到百分之百。
     我的算法是这样的，还是0-10长度为10的模拟。
     以0开始5结束（这也是我原先理解中的0-5），可以得到一条长度为5的msd-t曲线（也就是上面有5个点）。同理1-6,2-7,3-8,4-9,5-10，一共有6条msd-t曲线，每条长度都只有5，相当于30个点经过平均绘制出来了最终只有5个点的长度为5的msd曲线。
      对比一下，老师您说的算法，一共用到了55个点画出了一条长度为10的（有10个点）的曲线，确实数据利用率高多了！
      话说回来，我的算法倒是理论上每个点的精度都是一样的，应该不存在越往后曲线误差越大。
      再次感谢sob老师，我感觉这下应该理解更进了一步！