使vmd实时显示gromacs轨迹中原子的受力 | 宜武汇-ag真人国际厅网站

显示原子受力有时很重要,但实现起来又是个麻烦事。由于这两天要用这个功能,所以写了一个vmd的tcl脚本,可以用箭头显示指定原子的受力,并且随进度条拖动实时更新,在这里做一下介绍。

vmd当前版本并不读取轨迹中的力,不能在内部直接调用,所以需要手动把轨迹中的力提取出来并读入vmd。在此例中,我们要分析序号为977的原子的受力,序号由0开始算。

这里用到一个技巧,结合gmxdump和grep可以容易地从.trr文件中提取977号原子受力,虽然比较慢,但是很方便:

gmxdump -f all.trr |grep “f.  977” > force.txt

包含f[  977]的行就被提取出来了。force.txt内容如下:

应当注意,为了使每一行都对应一帧结构,.mdp的nstxout和nstfout应当相同。(最后一帧未必,因为最后一步的结构不管怎么设nstxout都会被输出,而这一帧未必会正好是输出力的那步,但这个问题无关紧要)

把force.txt放到tcl启动后的默认文件夹(比如d:\study\vmd186,在控制台输入pwd可显示),启动vmd,载入轨迹,在控制台运行:

我们已将全部受力读入一个数组f,例如f(x43)代表第43帧时977号原子受力的x分量。

然后运行下面的内容,定义一个名为showforce的过程,用于绘制箭头

运行:

trace variable vmd_frame(0) w showforce

说明跟踪vmd_frame(0)变量,这个变量每当id为0的分子的进度条移动后都会更新,内容是当前帧号。此时便触发showforce过程根据此帧的受力来重绘箭头。

最后定义要显示哪个原子受力,在此例中运行set atom 977

拖动进度条,就可以看到效果了,如图所示的红箭头,箭头的长度与受力大小成正比。

如果不想显示了,只需输入

trace vdelete vmd_frame(0) w showforce

graphics 0 delete all

如果又想显示,仍然是运行trace variable vmd_frame(0) w showforce

为了方便,我将上述内容写进了showforce.tcl脚本。简要来说使用方法如下:

1. 将force.txt和showforce.tcl放入vmd默认文件夹

2. 启动vmd,载入轨迹

3. 在控制台运行source showforce.tcl

4. set atom 977

即可。

showforce.tcl文件内容如下

ps:还可将此脚本改写,显示每一帧某原子的速度。

原文链接:http://sobereva.com/36

网络摘文,本文作者:15h,如若转载,请注明出处:https://www.15cov.cn/2023/08/27/使vmd实时显示gromacs轨迹中原子的受力/

发表评论

邮箱地址不会被公开。 必填项已用*标注

网站地图