计算化学公社

标题: 利用VMD进行轨迹分析的TCL脚本 [打印本页]
作者
Author: 丁越    时间: 2022-2-12 18:26
标题: 利用VMD进行轨迹分析的TCL脚本
本帖最后由 丁越 于 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]

  6. #calculate the geometric diatance of two selections
  7. set sum 0
  8. set var1 [atomselect top $sel1]
  9. set var2 [atomselect top $sel2]
  10. for {set i 0} {$i < $nframe} {incr i} {
  11.         $var1 frame $i
  12.         set val1 [measure center $var1] ;#[measure center $sel1 weight mass] can be used to calculate the center of mass
  13.         $var2 frame $i
  14.         set val2 [measure center $var2]
  15.         set dist [vecdist $val1 $val2]
  16.         set sum [expr $sum+$dist]
  17.         puts $myfile "$i\t $dist"
  18. }
  19. close $myfile
  20. puts "the average diatance is: [expr double($sum)/$nframe]"
  21. 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"

  7. plot \
  8.         '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（可以自行比对一下脚本差异，欢迎展示更多的写法）.
(, 下载次数 Times of downloads: 48)


  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. #=======================================

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

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

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

  49. #output the probability distribution
  50. foreach i [lsort -real -decreasing [array names num]] {
  51.         set p [expr double($num($i))/$nframe*100]
  52.         puts $myfile "$i\t $p"
  53. }        
  54. close $myfile
  55. 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. #**********************************

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

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

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

  35. #calculating the mean number of $sel along Z direction and output the calculated results to znumdist.dat
  36. foreach i [lsort -real -decreasing [array names num]] {
  37.         set Mnum($i) [expr double($num($i))/$nframe]
  38.         puts $myfile "$i\t $Mnum($i)"
  39. }
  40. 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. }
复制代码
统计的结果如下所示：
[attach]45117[/attach]

(, 下载次数 Times of downloads: 54)

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

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

(, 下载次数 Times of downloads: 50) (, 下载次数 Times of downloads: 53)














作者
Author: sobereva    时间: 2022-2-13 03:32
红色的shaded部分大抵是文中跑AIMD统计的水的氧的分布，只有液体的粒子分布曲线才会那么平滑
作者
Author: 丁越    时间: 2022-2-13 09:31
sobereva 发表于 2022-2-13 03:32
红色的shaded部分大抵是文中跑AIMD统计的水的氧的分布，只有液体的粒子分布曲线才会那么平滑

谢谢老师，我再跑跑轨迹了看结果怎么样
作者
Author: suihg    时间: 2022-2-13 10:02
学习了！！
作者
Author: DOUKUN    时间: 2022-3-1 20:51
本帖最后由 DOUKUN 于 2022-3-1 21:23 编辑

楼主太牛了，脚本很实用，请教沿着z方向统计某五元环的质心的数密度分布，因为是质心，不是一个实际的原子，不知道怎么选中合适的质心点，楼主有什么见解
作者
Author: 丁越    时间: 2022-3-1 22:42
DOUKUN 发表于 2022-3-1 20:51
楼主太牛了，脚本很实用，请教沿着z方向统计某五元环的质心的数密度分布，因为是质心，不是一个实际的原子 ...

measure center $sel1 weight mass返回值就是质心
作者
Author: DOUKUN    时间: 2022-3-2 08:24
丁越 发表于 2022-3-1 22:42
measure center $sel1 weight mass返回值就是质心

谢谢版主，我试一下
作者
Author: DOUKUN    时间: 2022-3-7 11:30
本帖最后由 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报错部分我用下划线标出来了，然后我对脚本的改动用红色字体标出来了

作者
Author: 丁越    时间: 2022-3-7 21:49
DOUKUN 发表于 2022-3-7 11:30
版主，我想测我体系中所有五元环的质心的数密度，这样套用了你的脚本，但是总会报错，您能看出问题来吗
se ...

我这两天太忙了，你先把http://bbs.keinsci.com/thread-157-1-1.html里面的基本语法看了就会写了。
一行一行的运行结果，看看哪里出错
作者
Author: DOUKUN    时间: 2022-3-8 07:58
丁越 发表于 2022-3-7 21:49
我这两天太忙了，你先把http://bbs.keinsci.com/thread-157-1-1.html里面的基本语法看了就会写了。
一行 ...

好的，谢谢，最近我也在看这个教程，因为编程基础薄弱，目前还没有看出问题来，你有时间了可以帮忙看看，实在忙不过来就算了，谢谢
作者
Author: 函数与激情    时间: 2023-2-8 21:01
楼主好，我最近用了下您计算原子密度的tcl脚本，发现好像有点问题，想请教下您。对整条轨迹帧数取平均的话每个bin内的原子数大概率是个小数，实际情况却非如此，使用您的脚本算出来一条动力学轨迹的每个bin内原子密度总是整数，似乎不太合理。检查了下代码，在第28行后加入puts "$zval($j)", 我发现zval($j)这个值是不变的呀(但我没看出来代码有问题)，所以我想最后的密度值似乎是最后一帧的原子密度而不是整个动力学轨迹的平均密度。附件是我测试的动力学轨迹，请教一下楼主我的测试是否有问题。



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3