博文
VMD加载GROMACS轨迹的速度信息与时间信息
|
2016-10-21 16:56:11
Sobereva曾给过一个tcl脚本, 用于解决VMD不能读入GROMACS轨迹速度信息的问题, 具体参考 使VMD读入Gromacs产生的trr轨迹中速度信息的方法. 最近我需要用到VMD的这个功能, 就稍微看了下tcl的语法, 在原代码的基础上改进了一点, 使其使用更简单. 同时结合Sobereva的另一段代码 使VMD播放轨迹时同步显示帧号, 在播放轨迹的同时显示出模拟时间.
下面模拟的是C20和C60分子的相撞过程, 播放轨迹时对每个原子根据速度大小进行着色.
运行MD前将grompp.mdp文件中trr与xtc的输出频率设为相同
使用gmx traj -f traj.trr -ov抽取traj.trr轨迹文件中的速度, 默认存为veloc.xvg
使用VMD加载初始的conf.gro文件和traj.xtc轨迹文件(直接使用trr文件可能更简单, 但速度稍慢). 也可直接使用命令vmd conf.gro traj.xtc
将vt.tcl脚本复制到轨迹文件所在目录下
VMD命令行窗口中执行source vt.tcl使脚本生效
VMD命令行窗口中执行loadveloc即可加载veloc.xvg文件中的速度. 如果速度文件的名称不是veloc.xvg, 则使用loadveloc 速度文件名即可
如果播放轨迹时需要显示时间, 可在VMD命令行窗口中执行showtime. 执行showtime off则关闭时间显示. 默认的时间间隔为0.5 fs, 起始时间为0, 如需更改, 可使用showtime on 时间间隔 起始时间
播放轨迹时对每个原子根据速度大小进行着色, 可通过Graphics | Representations... | Coloring Method | trajectory | Velocity. 如需根据某一方向速度大小着色, 可使用User(x方向), User2(y方向),User3(z方向)
更改颜色方案, 可使用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}{dtin0.5}{t0in0}}{global vmd_frame;global dt;global t0
set dt $dtinset t0 $t0inif{$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 [$Molnum]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}{$Molframe$iset Vx {};set Vy {};set Vz {}for{set j 0}{$j<$Natm}{incr j}{lappend Vx [lindex$txt[expr3*$j+1]]lappend Vy [lindex$txt[expr3*$j+2]]lappend Vz [lindex$txt[expr3*$j+3]]}$Molset vx $Vx$Molset vy $Vy$Molset vz $Vz$Molset user $Vx$Molset user2 $Vy$Molset user3 $Vzgets$Fxvg txt
set txt [split$txt\t]}close$Fxvg} |
https://blog.sciencenet.cn/blog-548663-1010146.html
上一篇：Martini实例教程：蛋白质
下一篇：陈旭盛: TCL培训教程
