|
|
本帖最后由 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 文件中"
|
|