计算化学公社

 找回密码 Forget password
 注册 Register
Views: 2727|回复 Reply: 5
打印 Print 上一主题 Last thread 下一主题 Next thread

[VMD] 五元环质心沿z轴数密度的tcl脚本问题

[复制链接 Copy URL]

22

帖子

0

威望

338

eV
积分
360

Level 3 能力者

跳转到指定楼层 Go to specific reply
楼主
之前参照一个帖子利用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

22

帖子

0

威望

338

eV
积分
360

Level 3 能力者

6#
 楼主 Author| 发表于 Post on 2022-3-14 21:49:57 | 只看该作者 Only view this author
本帖最后由 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 这个数组名字好像也不合适,因为这个时候数组的值还是一个坐标,感觉这个思路可能不好走通吧

22

帖子

0

威望

338

eV
积分
360

Level 3 能力者

5#
 楼主 Author| 发表于 Post on 2022-3-14 21:49:17 | 只看该作者 Only view this author
本帖最后由 DOUKUN 于 2022-3-15 08:16 编辑

谢谢帮助,谢谢指点

432

帖子

11

威望

3424

eV
积分
4076

Level 6 (一方通行)

4#
发表于 Post on 2022-3-14 18:41:50 | 只看该作者 Only view this author
DOUKUN 发表于 2022-3-14 14:53
谢谢谢谢,就是测量质心这个语句,因为它读出来是个坐标,不像atomselect top这种,读出来原子,下面就好 ...

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

22

帖子

0

威望

338

eV
积分
360

Level 3 能力者

3#
 楼主 Author| 发表于 Post on 2022-3-14 14:53:22 | 只看该作者 Only view this author
丁越 发表于 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]
还是不对,也不知道怎么改合适了,水平有限

432

帖子

11

威望

3424

eV
积分
4076

Level 6 (一方通行)

2#
发表于 Post on 2022-3-14 08:42:26 | 只看该作者 Only view this author
foreach j [[measure center $var  weight mass] list] {
                set zval($j) [[atomselect top "index $j"] get z]
这一段有问题,原本脚本的意思是foreach得到每帧原子序号的列表,然后由下面的一句get z得到z坐标值。写脚本时先把思路想好,不要生搬硬套
自由发挥,野蛮生长

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

GMT+8, 2024-11-24 13:40 , Processed in 0.197138 second(s), 22 queries , Gzip On.

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