VMD加载GROMACS轨迹的速度信息与时间信息

• 2016-10-21 16:56:11

Sobereva曾给过一个tcl脚本, 用于解决VMD不能读入GROMACS轨迹速度信息的问题, 具体参考 使VMD读入Gromacs产生的trr轨迹中速度信息的方法. 最近我需要用到VMD的这个功能, 就稍微看了下tcl的语法, 在原代码的基础上改进了一点, 使其使用更简单. 同时结合Sobereva的另一段代码 使VMD播放轨迹时同步显示帧号, 在播放轨迹的同时显示出模拟时间.

下面模拟的是C20和C60分子的相撞过程, 播放轨迹时对每个原子根据速度大小进行着色.

使用方法

1. 运行MD前将grompp.mdp文件中trr与xtc的输出频率设为相同
2. 使用gmx traj -f traj.trr -ov抽取traj.trr轨迹文件中的速度, 默认存为veloc.xvg
3. 使用VMD加载初始的conf.gro文件和traj.xtc轨迹文件(直接使用trr文件可能更简单, 但速度稍慢). 也可直接使用命令vmd conf.gro traj.xtc
4. 将vt.tcl脚本复制到轨迹文件所在目录下
5. VMD命令行窗口中执行source vt.tcl使脚本生效
6. VMD命令行窗口中执行loadveloc即可加载veloc.xvg文件中的速度. 如果速度文件的名称不是veloc.xvg, 则使用loadveloc 速度文件名即可
7. 如果播放轨迹时需要显示时间, 可在VMD命令行窗口中执行showtime. 执行showtime off则关闭时间显示. 默认的时间间隔为0.5 fs, 起始时间为0, 如需更改, 可使用showtime on 时间间隔 起始时间
8. 播放轨迹时对每个原子根据速度大小进行着色, 可通过Graphics | Representations... | Coloring Method | trajectory | Velocity. 如需根据某一方向速度大小着色, 可使用User(x方向), User2(y方向), User3(z方向)
9. 更改颜色方案, 可使用Graphics | Colors... | Color Scale | Method.

几点说明

• 这种基于tcl的方法可行, 但需要编写tcl脚本. 如果你不喜欢tcl脚本的话, 至少还有两个变通的替代方案: 1. 可以将gromacs的轨迹文件转换为lammps的轨迹文件, 因为vmd支持读取lammps轨迹文件中的速度. 2. 将速度写到pdb文件中的温度因子或电荷列中, 再根据相应的项进行着色. gmx traj的-cv选项一定程度上可以完成这点, 但效率太低.

• 对原子根据其速度进行着色并不总是最好的方法, 更好的方法是根据动能或温度进行着色. 一则原子类型不同时, 质量小的原子速度会相对较大, 二则人们对原子速度大小没有多少直观感觉, 添加颜色标尺时不易把握. 如果使用相对动能或温度, 就更加直观, 也更容易把握了. 只要对脚本稍加修改就可以达到这个目的.

• 对于温度或速度这种物理量而言, 使用分散颜色方案是最好的, VMD中接近这种颜色方案的是BWR, 但仍有不少差距. 如果想使用自定义的颜色方案, 除了自己写tcl代码以外, 目前我不知道有没有简单的办法. 有关发散颜色方案的信息可参考我的另外两篇博文: 数据展示：请选择好的颜色映射方案, 几种颜色映射方案的解析式.

• tcl中的trace可以使用回调函数(callback)进行变量跟踪, 使用trace时回调函数必须带有参数, 否则执行有问题. 此外, 新版本tcl中trace的语法有所改变, 上面脚本中的使用方式以后会废弃, 建议大家使用时尽量使用新版本的语法. 具体可参看下面的资料: VMD手册: Tcl callbacks, tcl手册: trace.

vt.tcl脚本

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

proc showtime { {opt on} {dtin 0.5} {t0in 0} } {
	global vmd_frame; global dt; global t0
	set dt $dtin
	set t0 $t0in
	if {$opt==on} {
		trace variable vmd_frame([molinfo top]) w traceframe
	} elseif {$opt==off} {
		draw delete all
		trace vdelete vmd_frame([molinfo top]) w traceframe
	}
}

proc traceframe {name elem ops} {
	global vmd_frame; global dt; global t0
 	draw delete all
	draw color white
	set time [format "%6.1f fs" [expr ($vmd_frame([molinfo top]) * $dt) + $t0]]
	draw text {0 -5 -5} "$time" size 4 thickness 4
}

proc loadveloc { {Fxvg veloc.xvg} } {
	set Mol  [atomselect top all]
	set Natm [$Mol num]
	set Nfrm [molinfo top get numframes]

	set Fxvg [open $Fxvg r]
	gets $Fxvg txt
	while {[string match {[\#@]*} $txt]} {
		gets $Fxvg txt
	}
	set txt [split $txt \t]

	for {set i 1} {$i<$Nfrm} {incr i} {
		$Mol frame $i
		set Vx {}; set Vy {}; set Vz {}
		for {set j 0} {$j<$Natm} {incr j} {
			lappend Vx [ lindex $txt [expr 3*$j+1] ]
			lappend Vy [ lindex $txt [expr 3*$j+2] ]
			lappend Vz [ lindex $txt [expr 3*$j+3] ]
		}

		$Mol set vx $Vx
		$Mol set vy $Vy
		$Mol set vz $Vz
		$Mol set user  $Vx
		$Mol set user2 $Vy
		$Mol set user3 $Vz

		gets $Fxvg txt
		set txt [split $txt \t]
	}

	close $Fxvg
}

评论

• 2016-10-22 10:21:14 李正 这两个分子相撞是怎么实现的啊？

• 2016-10-22 13:59:03 Jerkwin 把两个分子的速度设为方向相反, 运行NVE模拟即可.

• 2016-11-07 11:43:09 粥四文 gmx 不是内部命令？

• 2016-11-08 08:21:30 Jerkwin 按你安装好的gmx来写, 或许是gmx_mpi之类的

• 2016-11-08 11:06:11 粥四文 ”将vt.tcl脚本复制到轨迹文件所在目录下“，怎么弄？tcl文件，是下载一个tcl.pro软件编写吗？我是个初学者，不太理解。麻烦您了。
• 2016-11-08 21:46:11 Jerkwin 新建一个文本文档, 将下面的脚本复制粘贴进去, 保存为vt.tcl, 就可以了.