计算化学公社

 找回密码 Forget password
 注册 Register
Views: 6162|回复 Reply: 9

[VMD] 利用VMD进行轨迹分析的TCL脚本

[复制链接 Copy URL]

327

帖子

9

威望

1986

eV
积分
2493

Level 5 (御坂)

发表于 Post on 2022-2-12 18:26:53 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 丁越 于 2022-3-22 09:44 编辑

利用VMD进行轨迹分析的TCL脚本


  这些小脚本都是我目前边学边写的,由于不知道脚本考虑的周不周到,以及是否存在知识上的误区,遂贴出来让坛子里的各位老师把把关,解解惑

一、UNIX-like command

1. tcat
   该命令可以像linux下的cat一样展示文件内容
  1. proc tcat {filename} {
  2.         set fid [open $filename r]
  3.         while {[gets $fid line]>=0} {
  4.                 puts $line
  5.         }
  6.         close $fid
  7. }
复制代码
2. tgrep
   该命令可以像grep一样抓取文件内容
  1. proc tgrep { pattern filename} {
  2.         set f [open $filename r]
  3.         while {[gets $f line ]>=0} {
  4.                 if {[regexp $pattern $line]} {
  5.                         puts $line
  6.                 }
  7.         }
  8.         close $f
  9. }
复制代码
二、统计两个结构之间距离随时间的变化     
虽然该脚本ene老师在(http://bbs.keinsci.com/thread-14821-1-1.html)已经写过了,但是作为练习和发散思维的目的,用多种方式写出来也是挺好的。第一种方式和ene老师的类似,不过这里应用了更简洁的vecdist语法,可以直接返回两个矢量间的距离,类似的这样的语法还有vecsub、veclength等,可在VMD手册查看。该脚本的思路比较简单:定义两个selection,然后计算其几何中心之间的距离,写for循环循环所有帧即可(注意这里也可以定义质心,将measure center $var替换成measure center $sel weight mass即可)。最后可以附带个gnuplot绘图脚本,执行完计算距离后直接就将结果输出来了,如下所示。
  1. set sel1 "same fragment as {element H}"
  2. set sel2 "same fragment as {element S}"
  3. set nframe [molinfo top get numframes]
  4. set myfile [open Distance.dat w]

  5. #calculate the geometric diatance of two selections
  6. set sum 0
  7. set var1 [atomselect top $sel1]
  8. set var2 [atomselect top $sel2]
  9. for {set i 0} {$i < $nframe} {incr i} {
  10.         $var1 frame $i
  11.         set val1 [measure center $var1] ;#[measure center $sel1 weight mass] can be used to calculate the center of mass
  12.         $var2 frame $i
  13.         set val2 [measure center $var2]
  14.         set dist [vecdist $val1 $val2]
  15.         set sum [expr $sum+$dist]
  16.         puts $myfile "$i\t $dist"
  17. }
  18. close $myfile
  19. puts "the average diatance is: [expr double($sum)/$nframe]"
  20. exec gnuplot distance.gp
复制代码
  1. set grid
  2. set ylabel "Distance (A)"
  3. set xlabel "Frames"
  4. set term pngcario font "Times new roman, 8"
  5. set output "Distance.png"

  6. plot \
  7.         'Distance.dat' u 1:2 w l lw 2 lc rgb "red" t "Distance"
复制代码
第二种计算距离的方法可以用数组。已知两个点的坐标(x1,y1,z1)、(x2、y2、z2),那么距离就是|L|=sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2),所以应用数组计算两点距离的思路就是将x1、x2、y1...分别赋值给数组中的不同元素,然后利用距离公式求解就行。
  1. for {set i 0} {$i < $nframe} {incr i} {
  2.         $sel1 frame $i
  3.         set val1(x) [lindex [measure center $sel1] 0]
  4.         set val1(y) [lindex [measure center $sel1] 1]
  5.         set val1(z) [lindex [measure center $sel1] 2]
  6.         $sel2 frame $i
  7.         set val2(x) [lindex [measure center $sel2] 0]
  8.         set val2(y) [lindex [measure center $sel2] 1]
  9.         set val2(z) [lindex [measure center $sel2] 2]
  10.         
  11.         set dx2 [expr ($val1(x)-$val2(x))**2]
  12.         set dy2 [expr ($val1(y)-$val2(y))**2]
  13.         set dz2 [expr ($val1(z)-$val2(z))**2]
  14.         set dist [expr sqrt($dx2+$dy2+$dz2)]
  15.         
  16.         puts $myfile "$i\t $dist"
  17. }
复制代码
三、统计距离的分布
  该脚本的统计思路如下图所示,a、b是统计距离的最大最小值,tau为统计间隔。我们已经在上面的脚本中计算了两个selection的所有帧的距离,那么我们定义一个数组A然后将所有的距离赋值给数组的不同元素,紧接着我们循环数组A中的所有元素,判断如果元素i如果属于哪个距离间隔(如2<A(i)>3),那么我们就给统计元素个数的数组num加1,然后遍历循环所有的数组A元素就可以得出每个间隔tau内的数目分布。数组num中我们将元素定义为每个间隔的中点位置,值为统计的数目,取中点的目的是试想当间隔很小的时候,中值就可以近似等于该采样点的真实距离值,但实际情况下需要合理选择间隔大小,否则会造成数值问题。最后将统计值除以采样数目(帧数)就得到百分率了。类似的,角度分布、二面角分布也可以套用类似的思路去写。  ene老师的思路与该脚本有所不同,简单来说就是计算完某帧i的距离DIST后将第多少个tau间隔作为给统计数目数组的元素,然后其数组值加1(可以自行比对一下脚本差异,欢迎展示更多的写法).
num1.png


  1. #=======================================
  2. #                       Parameters
  3. #=======================================
  4. set sel1 "same fragment as {index 19}"
  5. set sel2 "same fragment as {index 2}"
  6. set nbin 100
  7. set myfile [open Dist_distribution.dat w]
  8. #=======================================

  9. #calculating the geometric diatance of two selections
  10. set nframe [molinfo top get numframes]
  11. set var1 [atomselect top $sel1]
  12. set var2 [atomselect top $sel2]
  13. for {set i 0} {$i < $nframe} {incr i} {
  14.         $var1 frame $i
  15.         $var1 update
  16.         set val1 [measure center $var1]
  17.         $var2 frame $i
  18.         $var2 update
  19.         set val2 [measure center $var2]
  20.         set dist($i) [vecdist $val1 $val2]
  21. }

  22. #searching for the Maxdist and Mindist
  23. set Mindist $dist(1)
  24. set Maxdist $dist(1)
  25. foreach i [array names dist] {
  26.         if {$dist($i)<$Mindist} {
  27.                 set Mindist $dist($i)
  28.         }
  29.         if {$dist($i)>$Maxdist} {
  30.                 set Maxdist $dist($i)
  31.         }
  32. }

  33. #calculating the number distribution
  34. set tau [expr ($Maxdist-$Mindist)/$nbin]
  35. for {set i 0} {$i<=$nbin} {incr i} {
  36.         set num([expr $Mindist+($i+0.5)*$tau]) 0
  37. }
  38. foreach i [array names dist] {
  39.         for {set j 0} {$j<$nbin} {incr j} {
  40.                 if {$dist($i)>=[expr $Mindist+$j*$tau] && $dist($i)<[expr $Mindist+($j+1)*$tau]} {
  41.                         set num([expr $Mindist+($j+0.5)*$tau]) [expr $num([expr $Mindist+($j+0.5)*$tau])+1]
  42.                 }
  43.         }
  44. }

  45. #output the probability distribution
  46. foreach i [lsort -real -decreasing [array names num]] {
  47.         set p [expr double($num($i))/$nframe*100]
  48.         puts $myfile "$i\t $p"
  49. }        
  50. close $myfile
  51. exec gnuplot Dist_distribution.gp
复制代码

四、统计沿着z方向原子数目分布  
该脚本是根据这篇文献写的(DOI: 10.1038/NMAT4793,文献我已经放到文末),目的是为了统计单位TiO2晶胞内沿着z方向上氧原子(H2O中的氧)的数目分布。简单说下代码的思路:基本思想和上面说过的一样,我们首先将z方向划分为nbin个间隔,接着我们通过get z得到每个氧原子的z坐标值,假如某个index的O原子的z坐标值落到了某个间隔内,那么我们给这个间隔统计数目的数组值加一,其数组元素为该区间的中值。然后,我们遍历所有帧中的氧原子最后再求其平均,就得到了z方向的平均数目分布。
  1. #**********************************
  2. #*                                                                    *                                 
  3. #*           Znumdist.tcl         *
  4. #*                                *
  5. #**********************************
  6. #revising parameters below manually!
  7. set sel "element O"
  8. set myfile [open znumdist.dat w]
  9. set nbin 200
  10. #**********************************

  11. set nframe [molinfo top get numframes]
  12. set size [measure minmax [atomselect top all]]
  13. set Zmin [lindex [lindex $size 0] 2]
  14. set Zmax [lindex [lindex $size 1] 2]
  15. set tau [expr ($Zmax-$Zmin)/$nbin]
  16. set var [atomselect top $sel]

  17. #calculating the sum of Z coordinate for each selected $sel
  18. for {set i 0} {$i < $nbin} {incr i} {
  19.         set num([expr $Zmin+($i+0.5)*$tau])  0
  20. }

  21. for {set i 0} {$i < $nframe} {incr i} {
  22.          $var frame $i
  23.          $var update
  24.         foreach j [$var list] {
  25.                 set zval($j) [[atomselect top "index $j"] get z]
  26.                 set k [expr floor([expr double(($zval($j)-$Zmin))/$tau])]
  27.                 set num([expr $Zmin+($k+0.5)*$tau]) [expr $num([expr $Zmin+($k+0.5)*$tau])+1]
  28.         }
  29.         puts -nonewline [format "      Calculating i = %d of %d ...\r" $i $nframe]
  30. }

  31. #calculating the mean number of $sel along Z direction and output the calculated results to znumdist.dat
  32. foreach i [lsort -real -decreasing [array names num]] {
  33.         set Mnum($i) [expr double($num($i))/$nframe]
  34.         puts $myfile "$i\t $Mnum($i)"
  35. }
  36. close $myfile
复制代码


在24~31行中求数目的代码中,还有另一种写法,就是通过一个区间一个区间去循环查询每个氧原子是否属于这个间隔内,但是存在的问题是计算量将是上面的N倍,不过作为练习的目的,尽可能尝试多种写法。
  1. for {set i 0} {$i < $nframe} {incr i} {
  2.          $var frame $i
  3.         foreach j [$var list] {
  4.                 set zval($j) [lindex [measure center [atomselect top "index $j"]] 2]
  5.                 for {set k 0} {$k <= $nbin} {incr k} {
  6.                         if {$zval($j) >= [expr $k*$tau+$Zmin] && $zval($j) < [expr ($k+1)*$tau+$Zmin]} {
  7.                         set num([expr $Zmin+($k+0.5)*$tau]) [expr $num([expr $Zmin+($k+0.5)*$tau])+1]
  8.                         }
  9.                 }
  10.         }        
  11. }
复制代码
统计的结果如下所示:

znum5000.png

  TiO2模型中最底层有8个O,再上一层有16个,这些都与上图中结果符合;距离TiO2表层越远的单位薄层内水分子个数平均为1,也与文献吻合。

但是这篇文献中我有两点不解之处还向各位老师请教:
1. 文中说模型包含1/4 ML subsurface Ti_int (1/4ML的次表层Ti的缺陷),模型如下所示。按这个说法,次表层一共16个Ti原子,难道是抠掉4个Ti原子么?感觉有点多啊,不知道我是不是理解错了。2. 文中给出的O原子分布是一个填色曲线,按理说z方向的数目应该是一个个的薄层内被统计原子的数目,是属于离散值,难道是将这些值进行展宽得到的么?

rutile.png 文献.png














Structure of a model TiO2 photocatalytic interface.pdf

1.79 MB, 下载次数 Times of downloads: 28

评分 Rate

参与人数
Participants 9
威望 +1 eV +38 收起 理由
Reason
1145075595 + 4 谢谢分享
DZW + 5
lonemen + 5 好物!
DOUKUN + 4 好物!
hdhxx123 + 5 好物!
卡开发发 + 5 谢谢分享
sobereva + 1
ZZU_SCU + 5 谢谢分享
ggdh + 5 谢谢分享

查看全部评分 View all ratings

自由发挥,野蛮生长

4万

帖子

99

威望

4万

eV
积分
89888

管理员

公社社长+计算化学玩家

发表于 Post on 2022-2-13 03:32:39 | 显示全部楼层 Show all
红色的shaded部分大抵是文中跑AIMD统计的水的氧的分布,只有液体的粒子分布曲线才会那么平滑
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办最高质量的各种计算化学类培训:初级量子化学培训班基础(中级)量子化学培训班分子动力学与GROMACS培训班量子化学波函数分析与Multiwfn程序培训班,内容介绍以及往届资料购买请点击链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的最佳途径。培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取培训最新消息、避免错过网上最有价值的计算化学文章!
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人,讨论范畴相同
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

327

帖子

9

威望

1986

eV
积分
2493

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2022-2-13 09:31:43 | 显示全部楼层 Show all
sobereva 发表于 2022-2-13 03:32
红色的shaded部分大抵是文中跑AIMD统计的水的氧的分布,只有液体的粒子分布曲线才会那么平滑

谢谢老师,我再跑跑轨迹了看结果怎么样
自由发挥,野蛮生长

21

帖子

0

威望

1416

eV
积分
1437

Level 4 (黑子)

发表于 Post on 2022-2-13 10:02:04 | 显示全部楼层 Show all
学习了!!

16

帖子

0

威望

128

eV
积分
144

Level 2 能力者

发表于 Post on 2022-3-1 20:51:43 | 显示全部楼层 Show all
本帖最后由 DOUKUN 于 2022-3-1 21:23 编辑

楼主太牛了,脚本很实用,请教沿着z方向统计某五元环的质心的数密度分布,因为是质心,不是一个实际的原子,不知道怎么选中合适的质心点,楼主有什么见解

327

帖子

9

威望

1986

eV
积分
2493

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2022-3-1 22:42:24 | 显示全部楼层 Show all
DOUKUN 发表于 2022-3-1 20:51
楼主太牛了,脚本很实用,请教沿着z方向统计某五元环的质心的数密度分布,因为是质心,不是一个实际的原子 ...

measure center $sel1 weight mass返回值就是质心
自由发挥,野蛮生长

16

帖子

0

威望

128

eV
积分
144

Level 2 能力者

发表于 Post on 2022-3-2 08:24:47 | 显示全部楼层 Show all
丁越 发表于 2022-3-1 22:42
measure center $sel1 weight mass返回值就是质心

谢谢版主,我试一下

16

帖子

0

威望

128

eV
积分
144

Level 2 能力者

发表于 Post on 2022-3-7 11:30:17 | 显示全部楼层 Show all
本帖最后由 DOUKUN 于 2022-3-7 11:33 编辑

版主,我想测我体系中所有五元环的质心的数密度,这样套用了你的脚本,但是总会报错,您能看出问题来吗
set sel1 "ringsize 5 from all"
set myfile [open znumdist.dat w]
set nbin 200
#**********************************

set nframe [molinfo top get numframes]
set size [measure minmax [atomselect top all]]
set Zmin [lindex [lindex $size 0] 2]
set Zmax [lindex [lindex $size 1] 2]
set tau [expr ($Zmax-$Zmin)/$nbin]
set var [measure center $sel1 weight mass]

#calculating the sum of Z coordinate for each selected $sel
for {set i 0} {$i < $nbin} {incr i} {
        set num([expr $Zmin+($i+0.5)*$tau])  0
}

for {set i 0} {$i < $nframe} {incr i} {
         $var frame $i
         $var update
        foreach j [$var list] {
                set zval($j) [[atomselect top "index $j"] get z]
                set k [expr floor([expr double(($zval($j)-$Zmin))/$tau])]
                set num([expr $Zmin+($k+0.5)*$tau]) [expr $num([expr $Zmin+($k+0.5)*$tau])+1]
        }
        puts -nonewline [format "      Calculating i = %d of %d ...\r" $i $nframe]
}

#calculating the mean number of $sel along Z direction and output the calculated results to znumdist.dat
foreach i [lsort -real -decreasing [array names num]] {
        set Mnum($i) [expr double($num($i))/$nframe]
        puts $myfile "$i\t $Mnum($i)"
}
set size [measure minmax [atomselect top all]]
set Zmin [lindex [lindex $size 0] 2]
set Zmax [lindex [lindex $size 1] 2]
set tau [expr ($Zmax-$Zmin)/$nbin]
set var [measure center $sel weight mass]

#calculating the sum of Z coordinate for each selected $sel
for {set i 0} {$i < $nbin} {incr i} {
        set num([expr $Zmin+($i+0.5)*$tau])  0
}

for {set i 0} {$i < $nframe} {incr i} {
         $var frame $i
         $var update
        foreach j [$var list] {
                set zval($j) [[atomselect top "index $j"] get z]
                set k [expr floor([expr double(($zval($j)-$Zmin))/$tau])]
                set num([expr $Zmin+($k+0.5)*$tau]) [expr $num([expr $Zmin+($k+0.5)*$tau])+1]
        }
        puts -nonewline [format "      Calculating i = %d of %d ...\r" $i $nframe]
}

#calculating the mean number of $sel along Z direction and output the calculated results to znumdist.dat
foreach i [lsort -real -decreasing [array names num]] {
        set Mnum($i) [expr double($num($i))/$nframe]
        puts $myfile "$i\t $Mnum($i)"
}
close $myfile
ringsize 5 from all
>Main< (crdensity) 16 % file31
>Main< (crdensity) 17 % 200
>Main< (crdensity) 18 % 1001
>Main< (crdensity) 19 % {0.0011849419679492712 0.0008726119995117188 6.0} {41.1982421875 41.999969482421875 146.260009765625}
>Main< (crdensity) 20 % 6.0
>Main< (crdensity) 21 % 146.260009765625
>Main< (crdensity) 22 % 0.701300048828125
>Main< (crdensity) 23 % measure center: no atom selection
>Main< (crdensity) 24 % >Main< (crdensity) 25 % can't read "var": no such variable
>Main< (crdensity) 26 % >Main< (crdensity) 27 % close $myfile报错部分我用下划线标出来了,然后我对脚本的改动用红色字体标出来了

327

帖子

9

威望

1986

eV
积分
2493

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2022-3-7 21:49:25 | 显示全部楼层 Show all
DOUKUN 发表于 2022-3-7 11:30
版主,我想测我体系中所有五元环的质心的数密度,这样套用了你的脚本,但是总会报错,您能看出问题来吗
se ...

我这两天太忙了,你先把http://bbs.keinsci.com/thread-157-1-1.html里面的基本语法看了就会写了。
一行一行的运行结果,看看哪里出错
自由发挥,野蛮生长

16

帖子

0

威望

128

eV
积分
144

Level 2 能力者

发表于 Post on 2022-3-8 07:58:03 | 显示全部楼层 Show all
丁越 发表于 2022-3-7 21:49
我这两天太忙了,你先把http://bbs.keinsci.com/thread-157-1-1.html里面的基本语法看了就会写了。
一行 ...

好的,谢谢,最近我也在看这个教程,因为编程基础薄弱,目前还没有看出问题来,你有时间了可以帮忙看看,实在忙不过来就算了,谢谢

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-2 22:20 , Processed in 0.608828 second(s), 26 queries .

快速回复 返回顶部 返回列表 Return to list