计算化学公社

标题: 五元环质心沿z轴数密度的tcl脚本问题 [打印本页]

作者
Author: DOUKUN 时间: 2022-3-13 18:04
标题: 五元环质心沿z轴数密度的tcl脚本问题
之前参照一个帖子利用VMD进行轨迹分析的TCL脚本 - 分子模拟 (Molecular Modeling) - 计算化学公社 (keinsci.com)内修改了一个脚本，大家帮忙看看哪里有问题，为什么还是算不出来，脚本哪里有问题呢？脚本内红字是我改动的地方
# measure center $sel1 weight mass
set sel "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 [atomselect top $sel]
#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 [[measure center $var  weight mass] 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

作者
Author: 丁越 时间: 2022-3-14 08:42
foreach j [[measure center $var weight mass] list] {
set zval($j) [[atomselect top "index $j"] get z]
这一段有问题，原本脚本的意思是foreach得到每帧原子序号的列表，然后由下面的一句get z得到z坐标值。写脚本时先把思路想好，不要生搬硬套

作者
Author: DOUKUN 时间: 2022-3-14 14:53

丁越发表于 2022-3-14 08:42
foreach j [[measure center $var weight mass] list] {
set zval($j) [[atomselect top ...

谢谢谢谢，就是测量质心这个语句，因为它读出来是个坐标，不像atomselect top这种，读出来原子，下面就好分析了。这一段确实有问题，但是也尝试过各种修改，放在别的位置，或者直接set COM [measure center $sel weight mass]
set zval($j) [lindex $COM 2]
还是不对，也不知道怎么改合适了，水平有限

作者
Author: 丁越 时间: 2022-3-14 18:41

DOUKUN 发表于 2022-3-14 14:53
谢谢谢谢，就是测量质心这个语句，因为它读出来是个坐标，不像atomselect top这种，读出来原子，下面就好 ...

你得想个法子让某个属性值代表整个五元环，比如这个五元环中仅含有一个氮原子，你可以sel选中当前帧中所有的N原子，然后建立个数组，数组元素为每个环上的N原子，数组值为五元环的质心，这样就可以统计分布了

作者
Author: DOUKUN 时间: 2022-3-14 21:49
本帖最后由 DOUKUN 于 2022-3-15 08:16 编辑

谢谢帮助，谢谢指点

作者
Author: DOUKUN 时间: 2022-3-14 21:49
本帖最后由 DOUKUN 于 2022-3-15 08:56 编辑

丁越发表于 2022-3-14 18:41
你得想个法子让某个属性值代表整个五元环，比如这个五元环中仅含有一个氮原子，你可以sel选中当前帧中所 ...

谢谢你的点播：
set zval($j) [[atomselect top "index $j"] get z]
set k [expr floor([expr double(($zval($j)-$Zmin))/$tau])]
实际上z坐标就是这一行所得，然后z就成为了数组zval，紧接着下面就运算了。z如何能改成质心的z呢？这部分多加一个数组，就是com（$sel）[measure center [atomselect top "ringsize 5 from index $j" weight mass]?这样好像也不对，需要去脚本最前面进行比较大的改动？建立好数组了，只确立了某个原子的坐标，转移到了其质心上面，但atomselect 这个数组名字好像也不合适，因为这个时候数组的值还是一个坐标，感觉这个思路可能不好走通吧

欢迎光临计算化学公社 (http://bbs.keinsci.com/)