均方位移的计算求助

北京的大猫

各位前辈及sob 老师，
      本人因为后续工作的需要，需自己计算msd（不使用 gmx msd）,对于gmx在周期性边界条件下得到的轨迹文件，使用trjconv命令，先-pbc res然后-pbc nojump得到脱去
PBC的轨迹文件，利用该轨迹文件得到的质心计算第0 ps到第1 ps的msd,发现与gmx msd（0-1ps）计算的有明显的差别(均未进行系综平均),请问msd应该怎么计算呢？


     谢谢各位老师的帮助！

sobereva

你先得说清楚你是怎么算的，别人才可能告诉你差异在哪
记得gmx msd会对轨迹中所有可以取的1ps跨度算的MSD做平均，不是只计算一次0~1 ps的MSD

北京的大猫

sobereva 发表于 2021-4-14 17:40
你先得说清楚你是怎么算的，别人才可能告诉你差异在哪
记得gmx msd会对轨迹中所有可以取的1ps跨度算的MSD ...

谢谢sob老师，
我是这样算的：

对于nojump得到的文件计算质心坐标，取0-1ps（绝对时间）的坐标，计算(sum [R(t=1)-R(0)]^2)/N，N为粒子数，与gmx的计算结果比较：
gmx  msd  -f   xxx.trr    -s   xxx.tpr  -b  0   -e  1  然后选择同种粒子，
但是 两者得到的结果明显不同，用同样的方法计算0-5ps,每秒的msd也是有明显不同，但又不是相差明显的倍数。

sobereva

北京的大猫 发表于 2021-4-14 19:31
谢谢sob老师，
我是这样算的：

把问题简化，先只计算一个原子看看是否相符，这样就没质心的事情

北京的大猫

sobereva 发表于 2021-4-15 07:40
把问题简化，先只计算一个原子看看是否相符，这样就没质心的事情

sob老师，对于一个原子的情况，自己写了一个索引组，gmx计算的结果和自己的结果是一致的，

对于一个含有多个原子的阴离子，对-res ,-nojump之后的PDB文件，计算某个离子0-1ps的质心，利用质心计算msd,
与利用gmx计算同一个阴离子0-1ps内的msd比较仍然不同，自己的质心程序验证过没有问题，难道说gmx计算多原子的阴离子的msd的不是按照质心来计算的吗?

lyj714

分子加-mol，默认单原子

北京的大猫

lyj714 发表于 2021-4-15 14:55
分子加-mol，默认单原子

谢谢您，用-mol，-nojump试过了，结果还是对不上。

lyj714

北京的大猫 发表于 2021-4-15 15:12
谢谢您，用-mol，-nojump试过了，结果还是对不上。

用gmx自带工具导出某一个分子的质心坐标变化，然后和你自己程序计算结果比对，我认为可能是这个原因，gmx自己对pbc的处理要全面，就算你导出的pdb/gro处理了轨迹，但是也可能你算的不对。
另外，我说的-mol结合的index必须也是分子编号哦，而不是默认的原子index

北京的大猫

lyj714 发表于 2021-4-15 15:20
用gmx自带工具导出某一个分子的质心坐标变化，然后和你自己程序计算结果比对，我认为可能是这个原因，gmx ...

谢谢您，但我对比的那一个离子一直在盒子内，不涉及边界条件的问题。我用gmx自带的命令试试。

北京的大猫

lyj714 发表于 2021-4-15 15:20
用gmx自带工具导出某一个分子的质心坐标变化，然后和你自己程序计算结果比对，我认为可能是这个原因，gmx ...

试过了，gmx traj给出的质心和我的程序计算的质心结果是一样的。

lyj714

北京的大猫 发表于 2021-4-15 15:43
试过了，gmx traj给出的质心和我的程序计算的质心结果是一样的。

那我倒是看不出还有啥原因，注意我刚刚的回复说了-mol结合的必须是分子index文件。我有时间晚上再写个程序测试下看看。还有个东西你必须注意，gmx msd的默认重启点时间间隔是10 ps，也就是说你0-1ps可能只有一帧，那可能结果本身就不对。

北京的大猫

lyj714 发表于 2021-4-15 15:45
那我倒是看不出还有啥原因，注意我刚刚的回复说了-mol结合的必须是分子index文件。我有时间晚上再写个程 ...

不明白-mol为什么要结合分子的index文件，gmx自己不就可以识别哪个分子吗？
另外，我算的就是0-1ps一个阴离子的一帧的msd，如果一帧的结果都对不上的话，后续的结果就更无法对了吧?（小白一枚，请多多指教）谢谢您。

lyj714

北京的大猫 发表于 2021-4-15 16:00
不明白-mol为什么要结合分子的index文件，gmx自己不就可以识别哪个分子吗？
另外，我算的就是0-1ps一个 ...

还有可能你参考也就是t=0的坐标是不是弄错了，gmx可是读的tpr中的，不知道你用的哪个。

北京的大猫

lyj714 发表于 2021-4-15 16:05
还有可能你参考也就是t=0的坐标是不是弄错了，gmx可是读的tpr中的，不知道你用的哪个。

即便gmx使用tpr中的坐标，在考虑浮点计算误差的前提下，两者的结果还是不同的。

lyj714

北京的大猫 发表于 2021-4-15 19:32
即便gmx使用tpr中的坐标，在考虑浮点计算误差的前提下，两者的结果还是不同的。

我简单测试了一下（下图为一个水分子），设置-trestart 0和自己写的代码计算结果基本一致。