DNA直径—RES2个P

注释：我现在需要一个vmd的tcl脚本进行分析直径为R随时间的变化，依次计算残基2内元素P之间的距离为D1，计算残基3内元素P之间的距离为D2，依次类推，计算残基12内元素P之间的距离为D11，直径R为D1，D2，....D11之和的平均数，依次将帧数T，D1，D2.... D11,R写入到radius.txt文件中，每列保留8个字符，其中小数点后5个字符，每列均为右对齐，每列之间空4个字符，其中第一行将T，D1，D2.... D11,R对应的写上。

# VMD脚本：计算核酸直径随时间变化

set residues {2 3 4 5 6 7 8 9 10 11 12}  ;# 分析残基范围

# 创建输出文件

set outfile [open "radius.txt" w]

# 写入标题行

set header "       T    "

for {set i 1} {$i <= 11} {incr i} {

    set header [format "%s      D%-1d    " $header $i]

}

set header [format "%s        R" $header]

puts $outfile $header

# 获取总帧数

set num_frames [molinfo top get numframes]

# 遍历每一帧

for {set frame 0} {$frame < $num_frames} {incr frame} {

    # 更新到当前帧

    animate goto $frame

    # 初始化变量

    set distances {}

    set total_sum 0.0

    # 遍历每个残基

    foreach res $residues {

        # 选择残基内的P原子

        set sel [atomselect top "resid $res and name P"]

        set atoms [$sel get index]

        set num_atoms [llength $atoms]

        # 计算残基内P原子最大距离

        set max_dist 0.0

        if {$num_atoms >= 2} {

            set coords [$sel get {x y z}]

            for {set i 0} {$i < $num_atoms} {incr i} {

                for {set j [expr {$i+1}]} {$j < $num_atoms} {incr j} {

                    set dist [vecdist [lindex $coords $i] [lindex $coords $j]]

                    if {$dist > $max_dist} {

                        set max_dist $dist

                    }

                }

            }

        }

        $sel delete

        lappend distances $max_dist

        set total_sum [expr {$total_sum + $max_dist}]

    }

    # 计算直径R（平均距离）

    set R [expr {$total_sum / [llength $residues]}]

      # 显示进度

            puts "已处理: [format %6d/$nframes $fr]帧..."

            flush stdout

    # 格式化输出行

    set line [format "%8d    " $frame]

    foreach dist $distances {

        set line [format "%s%8.5f    " $line $dist]

    }

    set line [format "%s%8.5f" $line $R]

    # 写入结果

    puts $outfile $line

}

# 关闭文件

close $outfile

puts "分析完成！结果已保存至 radius.txt"