24小时热门版块排行榜    

查看: 1942  |  回复: 15

jianchaoyv

金虫 (小有名气)

[交流] 【求助】向bay__gulf求助!! 已有4人参与

目前在做碳纳米管中水的相变,请教一下怎样写.tcl脚本来计算radial density of oxygen 和 orientation distributions of water molecules within CNTs?
虽然这个问题对你可能比较简单,但对我来说却很棘手,望赐教!
同时恭喜bay__gulf发JPCA!
也欢迎其他高手指教!
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

bay__gulf

金虫 (著名写手)

刘苏州

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
jianchaoyv(金币+10):谢谢指点 2010-06-22 18:29:11
ghcacj(金币+10):谢谢 2010-06-22 19:55:31
radial density, 也就是径向分布, 很常用, vmd 也提供了一个工具可以直接计算
就是measure gofr, 参考 http://www.ks.uiuc.edu/Research/vmd/current/ug/node133.html, 以及一个gui 界面, 参考 http://www.ks.uiuc.edu/Research/vmd/plugins/gofrgui/

至于取向分布, 用户较少只能自己写了,
vmd 的强大之处不在于内置了多少工具, 更不是做多漂亮的图形
而是提供了一个原子(集团)为数据结构的编程环境, 可以按照自己的需要实现任意的功能
下面是当时用的代码, 现在看来很幼稚了但拿它入门还是不错的

# count the dipolar-z angle and Legendre 2nd order polynomial in a pore,
# make sure the psf/dcd be the top molecular

set outf [open theta.dat w]

set numFrame [molinfo top get numframes]

set lgd 0
for { set fr 1 } { $fr < $numFrame } {incr fr } {
   set wt [atomselect top "water and name OH2 and z < 7.6 and z > -7.6 " frame $fr]
   set n [$wt num]
   set theta 0

   foreach resid [$wt get resid] {

      set wto  [atomselect top "resid $resid and name OH2" frame $fr]
      set wth1 [atomselect top "resid $resid and name H1" frame $fr]
      set wth2 [atomselect top "resid $resid and name H2" frame $fr]

      set coordo  [measure center  $wto]
      set coordh1 [measure center $wth1]
      set coordh2 [measure center $wth2]

      set coordh  [vecadd $coordh1 $coordh2]
      set coordh  [vecscale 0.5 $coordh]
      set vecd    [vecsub $coordh $coordo]

      set thetai  [expr acos ([expr [lindex $vecd 2] / [veclength $vecd]])]
      set thetai  [expr $thetai * 180.0 / 3.1416]
      set theta   [expr $theta + $thetai]

      $wto  delete
      $wth1 delete
      $wth2 delete
   }

   set theta [expr $theta/$n ]
   set costheta [expr cos ([expr $theta*3.1416/180.0])]
   set lgd [expr $lgd + [expr $costheta * $costheta]]
   puts $fr
   puts $outf [format "%d %f" $fr $theta]

   $wt delete
}
   set lgd [expr $lgd / $numFrame * 1.5 - 0.5]
   puts "3*<(cosθ)^2> - 1 = $lgd"
close $outf
2楼2010-06-22 16:44:02
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

jianchaoyv

金虫 (小有名气)

能否具体讲一下两种操作方法的步骤或给一个具体例子?你给的连接看了,但还是不知道怎么去操作。麻烦你帮忙写一个.tcl脚本,谢谢!

[ Last edited by jianchaoyv on 2010-6-22 at 20:24 ]
3楼2010-06-22 19:41:03
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jianchaoyv

金虫 (小有名气)

请问上面提供的脚本中:
  set coordh  [vecadd $coordh1 $coordh2]
      set coordh  [vecscale 0.5 $coordh]
      set vecd    [vecsub $coordh $coordo]

      set thetai  [expr acos ([expr [lindex $vecd 2] / [veclength $vecd]])]
还有最后set lgd [expr $lgd / $numFrame * 1.5 - 0.5]
这几句看不懂,能不能解释一下?谢谢!

[ Last edited by jianchaoyv on 2010-6-23 at 09:56 ]
4楼2010-06-23 09:50:22
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

bay__gulf

金虫 (著名写手)

刘苏州

★ ★ ★ ★ ★
ghcacj(金币+5):谢谢 2010-06-24 10:06:00
set coordh  [vecadd $coordh1 $coordh2] ;#两个H的坐标之和
set coordh  [vecscale 0.5 $coordh] ;#两个H的平均坐标, 即中点
set vecd    [vecsub $coordh $coordo] ;#O原子到H中点的矢量, 即水的偶记(方向)
set thetai  [expr acos ([expr [lindex $vecd 2] / [veclength $vecd]])]
# 水的偶极u跟z轴(0 0 1)的夹角, cos(theta) = u(z)/|u|

还有最后set lgd [expr $lgd / $numFrame * 1.5 - 0.5]
勒让德项 (3-1)/2,
参考 http://www.mdbbs.org/viewthread.php?tid=11530
5楼2010-06-24 09:13:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jianchaoyv

金虫 (小有名气)

请问vecadd、vecscale、vecsub、veclength这些命令在哪里可以学习一下?根本不知道它们的用途?请给一个关于这类命令的全面的学习资料!!!十分感谢!

另外,这句set thetai  [expr acos ([expr [lindex $vecd 2] / [veclength $vecd]])]的意思 #水的偶极u跟z轴(0 0 1)的夹角, cos(theta) = u(z)/|u|看不懂啊?
我的理解是:lindex $vecd 2 #是O原子的位置矢量
            veclength $vecd #O原子到H中点的矢量|u|
    所以,set thetai  [expr acos ([expr [lindex $vecd 2] / [veclength   $vecd]])] 并不是真正的#水的偶极u跟z轴(0 0 1)的夹角, cos(theta) = u(z)/|u|
正在研究你给的这个脚本中!

[ Last edited by jianchaoyv on 2010-6-25 at 16:12 ]
6楼2010-06-25 15:28:58
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

bay__gulf

金虫 (著名写手)

刘苏州

★ ★ ★
jianchaoyv(金币+3): 2010-06-26 08:49:39
lei0736(金币+3):谢谢 2010-06-26 09:17:42
请问vecadd、vecscale、vecsub、veclength这些命令在哪里可以学习一下?根本不知道它们的用途?请给一个关于这类命令的全面的学习资料!!!十分感谢!
====
vmd ug 的第10章讲的是Vectors and Matrices
这个ug 仅仅是vmd 特有的和常用的tcl 命令, 更深入了解需要Tcl 的资料

lindex $vecd 2 #是O原子的位置矢量
===
vecd 是偶极矢量,
lindex $vecd 2 是偶极的z 分量
veclength $vecd 是偶极的长度

很荣幸自己的工作被别人分析和引用
7楼2010-06-25 16:34:06
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

coolrainbow

木虫 (著名写手)

未来国家冻凉

★ ★
jianchaoyv(金币+2): 2010-06-26 08:49:58
lei0736(金币+2):谢谢 2010-06-26 09:18:03
引用回帖:
Originally posted by jianchaoyv at 2010-06-25 15:28:58:
请问vecadd、vecscale、vecsub、veclength这些命令在哪里可以学习一下?根本不知道它们的用途?请给一个关于这类命令的全面的学习资料!!!十分感谢!

另外,这句set thetai  [expr acos ([expr [lindex $vec ...

赞!bay_gulf 可是这里水平极高而又热心助人的高手之一啊

建议你可以用一天的时间学学TCL语言,再看看vmd中特有的命令,再学习几个namd中已有的tcl脚本作为范例,基本上你就可以自己写脚本了

纯tcl推荐学习:Practical Programming in Tcl and Tk, Fourth Edition
技术博客:http://hi.baidu.com/coolrainbow/blog
8楼2010-06-25 21:08:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jianchaoyv

金虫 (小有名气)

bay__gulf确实是牛人,也很耐心!
我把你那段程序的理解写下来:
程序开始到此已理解
set theta [expr $theta/$n ]  #对每帧所有水分子偶极距求平均,另外,这个theta角是跟Z轴正半轴的夹角吧?而看到的文献中大都是算的与Z轴负半轴的夹角,感觉没什么差别,这样理解多吗?
   set costheta [expr cos ([expr $theta*3.1416/180.0])]#求cos(theta)
   set lgd [expr $lgd + [expr $costheta * $costheta]]#用于求不同帧的cos(theta)^2求和
   puts $fr #输出第几帧
   puts $outf [format " %f"  $theta] 输出该帧下水分子偶极距的平均值
   $wt delete
}
   set lgd [expr $lgd / $numFrame * 1.5 - 0.5] #对所有帧的2阶勒让德函数求平均
   puts "3*<(cosθ)^2> - 1 = $lgd" #输出P2的平均值
close $outf

对上面的解释若有错误请指正,还想请教一下算P2用来表征什么的?不知道算P2的目的。
若set theta [expr $theta/$n ]  #对每帧所有水分子偶极距求平均,理解正确的话,也就默认了每帧里的所有水分子的偶极距方向都一样了,感觉同一帧不同位置的水分子偶极距还是有可能不同的,有这种可能吗?

[ Last edited by jianchaoyv on 2010-6-26 at 10:09 ]
9楼2010-06-26 09:48:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

bay__gulf

金虫 (著名写手)

刘苏州

★ ★ ★ ★
小木虫(金币+0.5):给个红包,谢谢回帖交流
lei0736(金币+3):谢谢 2010-06-26 11:43:25
引用回帖:
Originally posted by jianchaoyv at 2010-06-26 09:48:01:
bay__gulf确实是牛人,也很耐心!
我把你那段程序的理解写下来:
程序开始到此已理解
set theta [expr $theta/$n ]  #对每帧所有水分子偶极距求平均,另外,这个theta角是跟Z轴正半轴的夹角吧?而看到的文献中 ...

理解很正确
引用回帖:
也就默认了每帧里的所有水分子的偶极距方向都一样了,感觉同一帧不同位置的水分子偶极距还是有可能不同的,有这种可能吗?

为使程序好理解, 修改了点数据
water and name OH2 and z < 7.6 and z > -7.6
这是整个管子中的水
改变z 的上下限, 可以求得一小段区域中, 水的平均偶极方向
多做几个位置, 可以求得不同位置的数值

[ Last edited by bay__gulf on 2010-6-26 at 10:37 ]
10楼2010-06-26 10:36:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 jianchaoyv 的主题更新
普通表情 高级回复 (可上传附件)
信息提示
请填处理意见