vmd脚本求助 分析分子周围的粒子分布

diaok

这个分析粒子分布的脚本是仿照sob老师的氢键分析写的

1. proc rgwithin {sel sel2 fps1 fps2} {
   
2. set sum 0.0
   
3. set cut  6  
   
4. #  0.1:0.1:6.1
   
5. set piece  0.1
   
6. set in  0
   
7. for {set j 1} {$in<=$cut} {incr j} {
   
8. set in  [expr $j*0.1]
   
9. set selin [atomselect top  "$sel && within $in of $sel2"]
   
10. set sum  0
    
11. for {set i $fps1} {$i<=$fps2} {incr i} {
    
12. $selin frame $i
    
13. $selin update
    
14. set tmp [$selin  num]
    
15. set sum [expr $sum+$tmp]
    
16. }
    
17. #echo  $in
    
18. lappend rg    [expr $sum/[expr $fps2-$fps1+1.0]]
    
19. $selin  delete
    
20. }
    
21. return  $rg
    
22. }

复制代码

但是这个计算速度特别慢，还会导致电脑特别卡
我的体系里面分析600帧就用了6分钟
而且不知道是什么原因进行第二次计算会越来越慢
脚本中明明已经用了$selin  delete

用什么方法计算一个原子到另一个大分子的最近距离比较快？
用gromacs或者vmd来实现都行

ruanyang

tcl脚本计算本省就比较慢，你改成向量计算可能会快一点！

sobereva

把轨迹转换成xyz格式，自行用Fortran或C代码来载入和分析会很快
注意运行过程中看看内存占用率，有没有越来越高的问题

diaok

sobereva 发表于 2015-12-4 18:23
把轨迹转换成xyz格式，自行用Fortran或C代码来载入和分析会很快
注意运行过程中看看内存占用率，有没有越 ...

sob老师回复神速啊

第二次提取的时候主要是cpu占用升到100%了
vmd的内存占用100m一直没怎么变

很奇怪批处理为什么会越来越慢，6分钟一次还算可以接受
而且删除了分子重新载入都没效果，只能重启vmd


一次批处理的代码

1. set outfile [open {W:\susedata\kbff\na8.txt} w]
   
2. set  tmp  [rgwithin   "type NA"  "resname BGC"   0 499  ]
   
3. puts $outfile $tmp
   
4. set  tmp  [rgwithin   "type OO"  "resname BGC"   0 499   ]
   
5. puts $outfile $tmp
   
6. set  tmp  [rgwithin   "type HH"  "resname BGC"   0 499  ]
   
7. puts $outfile $tmp
   
8. set  tmp  [rgwithin   "type OW"  "resname BGC"   0 299     ]
   
9. puts $outfile $tmp
   
10. close $outfile

复制代码

Fortran和c语言不怎么熟悉
用matlab的textscan可行吗？轨迹500m，有30000个frame

diaok

ruanyang 发表于 2015-12-4 18:16
tcl脚本计算本省就比较慢，你改成向量计算可能会快一点！

有什么好的方法可以在gromacs或者vmd里实现吗？
gromacs我只想到求每个原子与大分子的所有距离，求最小值
然后对所有帧所有原子做直方图

sobereva

gromacs里自带了g_mindist工具，可以看看是否合用

diaok

sobereva 发表于 2015-12-4 22:36
gromacs里自带了g_mindist工具，可以看看是否合用

非常感谢

批量处理可以提取出来

1. mkdir oho
   
2. file=oho/oho.ndx
   
3. cat index.ndx >  $file
   
4. for i in `seq 1 45`
   
5. do
   
6. echo [oho$i] >> $file
   
7. ((t=$i+4672-2))
   
8. echo $t >> $file
   
9. done
   
10. #批量生成所有离子的group
    
11. 
    
12. for i in `seq 1 45`
    
13. do
    
14. gmx mindist -s md.tpr  -f md.xtc  -n oho/oho.ndx -od oho/o$i.xvg  << EOF
    
15. oho$i
    
16. 2   #用字母匹配会出错，原因未知
    
17. EOF
    
18. done
    
19. # 批量提取每个粒子到大分子的最小距离
    
20. 
    
21. 获得了所有距离
    
22. 剩下的就交给matlab去循环读取了

复制代码

ruanyang

强，脚本写的很牛逼

compume

好厉害，我是刚入行菜鸟一枚，想问一下，直接用python之类的语言我应该怎么从trr中读取出每个分子或者原子的位置信息呢？
以为trr是二进制文件，不知道怎么看啊……

sobereva

compume 发表于 2015-12-17 21:40
好厉害，我是刚入行菜鸟一枚，想问一下，直接用python之类的语言我应该怎么从trr中读取出每个分子或者原子 ...

自己用python读trr很难实现，建议你把trr用VMD转换为xyz轨迹，是二进制文件，这种轨迹格式非常简单，之后自己写程序分析处理很方便。当然体积比trr大不少。

compume

sobereva 发表于 2015-12-18 00:05
自己用python读trr很难实现，建议你把trr用VMD转换为xyz轨迹，是二进制文件，这种轨迹格式非常简单，之后 ...

THX，可是我的trr文件有几个G……这样处理岂不会很慢

sobereva

compume 发表于 2015-12-18 10:01
THX，可是我的trr文件有几个G……这样处理岂不会很慢

没法避免

compume

sobereva 发表于 2015-12-18 12:18
没法避免

好嘞，多谢sob老师~

zhangj

您好，请问sob老师氢键分析的脚本在哪里可以找到，非常感谢

tjuchan

没有考虑周期边界是不是会少很多氢键？