利用gmx读取LAMMPS模拟的轨迹并进行后处理分析

wangyueda · 发表于 Post on 2024-5-30 22:41:00

本帖最后由 wangyueda 于 2024-6-3 23:05 编辑

最近研究了一下如何利用gmx读取LAMMPS（简称lmp）模拟的轨迹并进行后处理分析，具体总结如下：

前期工作：

0.1 跑lmp的轨迹格式：其他版本的lmp我不太清楚，我这里使用的2023 Aug 2是可以保存成gmx可读的xtc格式的，具体in文件中相关参数如下：

dump 1 all xtc 100 dump.xtc

复制代码

这里我设置每隔100步就保存一次轨迹文件。

1. 相关软件版本

LAMMPS 2023 Aug 2
GROMACS 2021.7 (ps 这里我试了2022.6版本貌似有bug，具体我记得公社有个帖子说了这个问题，反正gmx>=2022貌似都有问题，哪位老师同学记得可以放在评论区)
OVITO 3.10.6
Multiwfn 3.8(dev), release date: 2024-May-8

2. gro文件和ndx文件获取

2.1 先用ovito将lmp的初始结构文件（一般为.data或者.lmp）转换为xyz文件，这里假设为init.xyz。再利用Multiwfn将init.xyz转换为gro文件init.gro，
具体命令为：

Multiwfn init.xyz
\n
100
\n
2
\n
34
\n
\n
ctrl c

复制代码

其中\n表示回车，这样就得到了init.gro，检查init.gro最后一行盒子尺寸是否正确。

2.2 利用

gmx make_ndx -f init.gro -o system.ndx

复制代码

命令制作体系的索引文件，一般使用a+元素名称即可（e.g. a C表示选择体系中的碳原子），选择完毕后输入q保存。（熟悉gmx的同学这步应该不需要我教）

3. 计算msd和rdf

3.1 计算msd，使用

gmx msd -s init.gro -f dump.xtc -n system.ndx -rmcomm -o msd.xvg

复制代码

计算msd，根据提示选择目标原子，这里如果想计算某一方向上的msd也可以，例如前面命令中加入-type x即为计算x方向的msd，计算完会在msd.xvg的第20行得到扩散系数。如果需要指定拟合区间可以使用-beginfit和-endfit指定。-rmcomm一定得加上，lmp貌似没有移除质心位移。

3.2 计算rdf, 使用

gmx rdf -s init.gro -f dump.xtc -n system.ndx -rmax 1

复制代码

根据提示选择计算原子和参考原子，还可以加入-cn顺带计算配位数，这里不再赘述

第一次写教程，本人水平有限，请各位老师多批评指正~

xiaowu759 · 发表于 Post on 2024-5-31 06:46:35

谢谢分享！

wangyueda · 发表于 Post on 2024-5-31 16:45:02

其实这里我有个小问题想请教下卢老师@sobereva ，就是这里gmx msd -rmcomm选项是从哪个文件中提取体系原子的质量信息的（没有质量的话自然没法算质心的位置）？我查了gmx手册没有查到，而这里的xtc和gro文件显然是不含有质量信息的...

sobereva · 发表于 Post on 2024-6-1 10:06:47

wangyueda 发表于 2024-5-31 16:45
其实这里我有个小问题想请教下卢老师@sobereva ，就是这里gmx msd -rmcomm选项是从哪个文件中提取体系原子 ...

可能是根据原子名猜的

wangyueda · 发表于 Post on 2024-6-3 23:00:42

本帖最后由 wangyueda 于 2024-6-3 23:02 编辑

sobereva 发表于 2024-6-1 10:06
可能是根据原子名猜的

好的，谢谢卢老师，我之前担心这么操作会没法读取元素质量信息从而导致-rmcomm“形同虚设”，后面特地用sobtop生成了一个体系的伪top文件（并不是真正用来跑gmx mdrun，只是单纯的想把元素质量信息给出），结合一个mdp文件使用gmx grompp创建出体系的tpr文件，想着生成的tpr文件肯定包含体系中各个元素的质量信息，结果用

gmx msd -s init.gro -f dump.xtc -n system.ndx -rmcomm -o msd.xvg

复制代码

和

gmx msd -s init.tpr -f dump.xtc -n system.ndx -rmcomm -o msd.xvg

复制代码

计算得到的扩散系数完全一样（包括拟合误差），这才觉得第一条指令gmx确实是读取了（或者是猜了）原子质量的（ps 同时了对比了加了-rmcomm和不加-rmcomm的结果，发现得到的扩散系数不一致）

zbz-ecl · 发表于 Post on 2024-6-16 17:14:13

楼主您好，我用您的方法生成的rdf图有点奇怪呈现锯齿状，请问您生成的rdf是正常的吗，希望能得到解答。附上我用的.gro以及lammps输出文件与输入结构

solution.gro (12.49 KB, 下载次数 Times of downloads: 3)

input.lammps (1.01 KB, 下载次数 Times of downloads: 2)

Gonglinquan · 发表于 Post on 2025-1-23 23:54:07

您好，看了您这篇帖子于是尝试着想去使用gmx计算msd，我用vmd读取data文件并输出gro文件，使用dump输出的xtc文件，我也有一个关于masses的疑惑，就是在lammps中一般是使用数字去表示的原子类型，其输出的gro文件在原子类型一列也确实是数字，这导致我在使用gmx计算msd时，出现报错“Masses were requested”，所以您通过ovito将data文件转换为xyz时，再通过multiwfn将其转换为gro文件时，这一步骤和vmd直接得到的gro文件有什么不同吗？

Lacrimosa · 发表于 Post on 2025-1-24 10:40:41

Gonglinquan 发表于 2025-1-23 23:54
您好，看了您这篇帖子于是尝试着想去使用gmx计算msd，我用vmd读取data文件并输出gro文件，使用dump输出的x ...

你在VMD里直接用下面的命令转换一下原子名就行了
[atomselect top "name 1"] set name C

goldNAN · 发表于 Post on 2025-1-24 17:10:54

wangyueda 发表于 2024-5-31 16:45
其实这里我有个小问题想请教下卢老师@sobereva ，就是这里gmx msd -rmcomm选项是从哪个文件中提取体系原子 ...

考虑用我这个脚本。https://github.com/tamaswells/La ... /master/data2gmx.py

Gonglinquan · 发表于 Post on 2025-1-25 04:19:25

Lacrimosa 发表于 2025-1-24 10:40
你在VMD里直接用下面的命令转换一下原子名就行了
[atomselect top "name 1"] set name C

谢谢您的回复，我去试试看

coookies · 发表于 Post on 2025-3-21 11:00:32

您好，请问计算msd的方法中-rmcomm是不是现在不能用了？在官网manual搜gmx msd没有关于-rmcomm的说明，直接使用时也发生了”Error in user input:Invalid command-line options Unknown command-line option -rmcomm“报错。不带这个option会不会导致计算不准确？

chenym · 发表于 Post on 2025-3-27 11:04:20

“如何利用gmx读取LAMMPS（简称lmp）模拟的轨迹并进行后处理分析”，可是为什么过程看起来像是用初始结构进行格式转化后重新计算得到相关计算结果，并不是直接处理LAMMPS得到的轨迹文件呢（我不太懂GROMACS，但是想找一个可以处理LAMMPS输出轨迹文件melt.xyz的软件计算MSD）

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[GROMACS] 利用gmx读取LAMMPS模拟的轨迹并进行后处理分析

评分 Rate

浏览过的版块