计算化学公社

标题: 均方位移的计算求助 [打印本页]

作者
Author: 北京的大猫 时间: 2021-4-14 17:23
标题: 均方位移的计算求助
本帖最后由北京的大猫于 2021-4-14 17:25 编辑

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

谢谢各位老师的帮助！

作者
Author: sobereva 时间: 2021-4-14 17:40
你先得说清楚你是怎么算的，别人才可能告诉你差异在哪
记得gmx msd会对轨迹中所有可以取的1ps跨度算的MSD做平均，不是只计算一次0~1 ps的MSD

作者
Author: 北京的大猫 时间: 2021-4-14 19:31

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也是有明显不同，但又不是相差明显的倍数。

作者
Author: sobereva 时间: 2021-4-15 07:40

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

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

作者
Author: 北京的大猫 时间: 2021-4-15 14:31
本帖最后由北京的大猫于 2021-4-15 14:38 编辑

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

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

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

作者
Author: lyj714 时间: 2021-4-15 14:55
分子加-mol，默认单原子

作者
Author: 北京的大猫 时间: 2021-4-15 15:12

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

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

作者
Author: lyj714 时间: 2021-4-15 15:20
本帖最后由 lyj714 于 2021-4-15 15:37 编辑

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

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

作者
Author: 北京的大猫 时间: 2021-4-15 15:27
本帖最后由北京的大猫于 2021-4-15 15:32 编辑

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

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

作者
Author: 北京的大猫 时间: 2021-4-15 15:43

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

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

作者
Author: lyj714 时间: 2021-4-15 15:45
本帖最后由 lyj714 于 2021-4-15 15:50 编辑

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

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

作者
Author: 北京的大猫 时间: 2021-4-15 16:00

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

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

作者
Author: lyj714 时间: 2021-4-15 16:05

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

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

作者
Author: 北京的大猫 时间: 2021-4-15 19:32

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

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

作者
Author: lyj714 时间: 2021-4-15 22:14
本帖最后由 lyj714 于 2021-4-15 22:20 编辑

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

我简单测试了一下（下图为一个水分子），设置-trestart 0和自己写的代码计算结果基本一致。 (, 下载次数 Times of downloads: 89)

作者
Author: 北京的大猫 时间: 2021-4-16 10:21
本帖最后由北京的大猫于 2021-4-16 10:22 编辑

lyj714 发表于 2021-4-15 22:14
我简单测试了一下（下图为一个水分子），设置-trestart 0和自己写的代码计算结果基本一致。

谢谢，gmx 计算的不是质心的MSD, 对于一个包含多个原子的group,会计算组内每个原子的位移，然后按原子量为权重计算group的(R(t)-R(0)^2 )。

作者
Author: lyj714 时间: 2021-4-16 10:24
本帖最后由 lyj714 于 2021-4-16 10:29 编辑

北京的大猫发表于 2021-4-16 10:21
谢谢，gmx 计算的不是质心的MSD, 对于一个包含多个原子的group,会计算组内每个原子的位移，然后按原子量 ...

我说了，-mol就是计算的分子质心然后算的质心msd，不用-mol就是当成单个原子对待。当然默认肯定都是带原子质量权重的

作者
Author: 北京的大猫 时间: 2021-4-16 10:57

lyj714 发表于 2021-4-16 10:24
我说了，-mol就是计算的分子质心然后算的质心msd，不用-mol就是当成单个原子对待。当然默认肯定都是带原 ...

恩恩.谢谢

作者
Author: 五月雨 时间: 2022-4-19 20:30

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

请问楼主在计算时是如何区分体系中不同的分子的？

作者
Author: sobereva 时间: 2022-4-27 07:16

五月雨发表于 2022-4-19 20:30
请问楼主在计算时是如何区分体系中不同的分子的？

建立索引文件，把关注的分子定义成一个组，计算时选择恰当的组

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