DNA直径-相反匹配

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

# 选择主分子

set mol [molinfo top]

# 获取总帧数

set nframes [molinfo $mol get numframes]

# 打开输出文件

set outfile [open "radius.txt" w]

# 写入标题行

set header [format "%8s" "T"]

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

    append header "    [format "%8s" "D$i"]"

}

append header "    [format "%8s" "R"]"

puts $outfile $header

# 处理每一帧

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

    # 更新到当前帧

    animate goto $frame

    # 初始化变量

    set distances {}

    set sum 0.0

    # 计算D1到D11

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

        # 计算残基索引

        set resA [expr {2 + $i}]

        set resB [expr {24 - $i}]

        # 选择磷原子

        set selA [atomselect top "resid $resA and name P"]

        set selB [atomselect top "resid $resB and name P"]

        # 计算距离

        if {[$selA num] == 1 && [$selB num] == 1} {

            set dist [vecdist [lindex [$selA get {x y z}] 0] \

                             [lindex [$selB get {x y z}] 0]]

        } else {

            set dist 0.0

            puts "Warning: Missing P atom in frame $frame at residue pair A$resA/B$resB"

        }

        # 存储结果

        lappend distances $dist

        set sum [expr {$sum + $dist}]

        # 清理选择

        $selA delete

        $selB delete

    }

    # 计算平均值R

    set R [expr {$sum / 11.0}]

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

            flush stdout

    # 构建输出行（确保每列8字符，保留5位小数）

    set line [format "%8d" $frame]  ;# 帧数用整数格式

    foreach dist $distances {

        append line "    [format "%8.5f" $dist]"

    }

    append line "    [format "%8.5f" $R]"

    # 写入文件

    puts $outfile $line

}

# 关闭文件

close $outfile

puts "Analysis complete. Results saved to radius.txt"