计算化学公社

标题: VMD的TCL脚本运行中报错求助 [打印本页]

作者
Author: nobody 时间: 2023-9-7 23:39
标题: VMD的TCL脚本运行中报错求助
本帖最后由 nobody 于 2023-9-8 10:19 编辑

各位老师好，我想用TCL语言写一个VMD程序脚本，计算Gromacs产生的xtc文件中每个TBP残基上的C1和C10原子组成的向量，求所获得的向量与z方向的夹角沿着z轴的数密度分布。由于实在不会写脚本，因此用AI生成了一个。在运行过程中报错cannot find attribute '1': in atomsel: get: ，请各位老师指点一下

# 创建一个列表来存储每个TBP残基上的几何中心连线向量
set centerline_vectors_list {}

# 创建一个列表来存储每一帧的夹角
set angle_list {}

# 遍历每一帧
for {set frame 0} {$frame < [molinfo top get numframes]} {incr frame} {
# 设置当前帧
animate goto $frame

# 获取每个TBP残基上的C1和C10原子的坐标
set sel [atomselect top "resname TBP and (name C1 or name C10)"]
set num_atoms [$sel num]
set frame_vectors {}

for {set i 1} {$i < $num_atoms} {incr i} {
      set atom_coord [$sel get $i]
      lappend frame_vectors $atom_coord
}

# 计算每个TBP残基的几何中心连线向量
if {[llength $frame_vectors] == 2} {
      set c1_coord [lindex $frame_vectors 0]
      set c2_coord [lindex $frame_vectors 1]
      set center_coord [vecscale 0.5 [vecadd $c1_coord $c2_coord]]
      set centerline_vector [vecsub $c2_coord $c1_coord]

      # 计算向量与z轴的夹角，并将角度保存到列表中
      set z_vector {0.0 0.0 1.0}
      set cos_angle [vecdot $centerline_vector $z_vector]
      set angle [expr acos($cos_angle) * 180.0 / 3.1415926]
      lappend angle_list $angle
}

# 关闭选择
$sel delete
}

# 绘制夹角数密度分布
set num_bins 180
set bin_width [expr 180.0 / $num_bins]
set hist [array create hist]
foreach angle $angle_list {
set bin [expr int($angle / $bin_width)]
set hist($bin) [expr $hist($bin) + 1]
}

# 输出数据到文件
set output_file [open "angle_distribution.dat" w]
foreach bin [array names hist] {
set angle_min [expr $bin * $bin_width]
set angle_max [expr ($bin + 1) * $bin_width]
set count $hist($bin)
puts $output_file "$angle_min-$angle_max $count"
}
close $output_file

puts "夹角数密度分布已保存到 angle_distribution.dat 文件中"

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