24小时热门版块排行榜    

查看: 2449  |  回复: 15
当前只显示满足指定条件的回帖,点击这里查看本话题的所有回帖

jianchaoyv

金虫 (小有名气)

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

目前在做碳纳米管中水的相变,请教一下怎样写.tcl脚本来计算radial density of oxygen 和 orientation distributions of water molecules within CNTs?
虽然这个问题对你可能比较简单,但对我来说却很棘手,望赐教!
同时恭喜bay__gulf发JPCA!
也欢迎其他高手指教!
回复此楼
已阅   回复此楼   关注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的回帖
查看全部 16 个回答

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的回帖
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考博] 售SCI一区T0P文章,我:8.O55.1.O.54,科目全,可十急 +3 2JOx3r2CYEgw 2026-08-21 8/400 2026-08-23 00:28 by OEbVnUOu01ol
[论文投稿] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 QTy3jDtz1uLt 2026-08-21 8/400 2026-08-22 23:51 by OEbVnUOu01ol
[硕博家园] 售SCI一区T0P文章,我:8O.55.1.O.54,科目全,可伽急 +4 QTy3jDtz1uLt 2026-08-21 10/500 2026-08-22 23:42 by OEbVnUOu01ol
[基金申请] 2026国自然函评费到账 +11 羊腰板 2026-08-21 11/550 2026-08-22 22:03 by lsbin3733
[基金申请] 93BebMhtakh前后11位开头都是大写 +7 且听虎啸 2026-08-17 8/400 2026-08-22 21:55 by 医学老男孩
[硕博家园] 售SCI一区T0P文章,我:8.O.55.1.O.5.4,科目全,可+急 +3 2JOx3r2CYEgw 2026-08-22 5/250 2026-08-22 21:54 by OEbVnUOu01ol
[教师之家] 跳槽后在研项目怎么办? +3 简单化xn 2026-08-22 7/350 2026-08-22 19:46 by 简单化xn
[基金申请] 什么时候开奖? +7 CrisMessi 2026-08-18 8/400 2026-08-22 18:29 by 淀粉搬运工
[基金申请] 只有每年这种时候来逛逛小木虫 +24 yaoyewhu2008 2026-08-20 26/1300 2026-08-22 17:43 by kammury
[基金申请] filecode,4个jtjc了 +13 ziyangfang 2026-08-19 16/800 2026-08-22 17:08 by WH3796
[教师之家] 售一区SCI文章T0P,我:8O.551.O54,科目全,可十急 +4 2JOx3r2CYEgw 2026-08-21 4/200 2026-08-22 16:52 by sunzitan
[基金申请] 人气不行了 +8 fansofjerry 2026-08-21 8/400 2026-08-22 16:30 by zyqchem
[基金申请] 今天基金会出结果吗?20260819 +16 kkkl_v 2026-08-19 17/850 2026-08-22 16:12 by 阿布Abu
[基金申请] 放榜前的不淡定 40+4 snowwithsea 2026-08-19 14/700 2026-08-21 23:51 by cratir
[基金申请] 建议基金发布提前给出明确的时间点 +10 kulium 2026-08-21 13/650 2026-08-21 21:58 by alongwaytogo
[基金申请] 感觉是下周放榜了 +7 angus9576 2026-08-17 12/600 2026-08-21 13:38 by weiyin
[基金申请] 基金啊基金 +4 longfie172 2026-08-20 4/200 2026-08-21 08:58 by mark mao
[基金申请] 时间戳变了,能看出什么问题? +18 基诺咪客 2026-08-17 23/1150 2026-08-20 17:19 by Godzela
[基金申请] 朋友圈看到的 +6 wangzilk 2026-08-18 8/400 2026-08-19 10:55 by Haru815
[基金申请] 今天维护系统维护 祝所有人 高中 +8 gjjjzhong 2026-08-18 9/450 2026-08-18 13:01 by 家与远方
信息提示
请填处理意见