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对应的写上。

# 计算DNA直径随时间变化的VMD Tcl脚本 - 精确配对版本

# 参数设置

set outfile "radius.txt" ;# 输出文件名

# 主程序

proc calc_dna_radius {{mol top}} {

    global outfile

    # 打开输出文件

    set fout [open $outfile w]

    # 写入标题行

    puts -nonewline $fout [format "%8s" "T"]

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

        puts -nonewline $fout [format "%12s" "D$i"]

    }

    puts $fout [format "%12s" "R"]

    # 获取总帧数

    set nframes [molinfo $mol get numframes]

    puts "开始分析DNA直径变化: 共 [format %d $nframes] 帧"

    # 遍历每一帧

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

        # 输出当前帧号

        puts -nonewline $fout [format "%8d" $fr]

        # 初始化变量

        set distances {}

        set sum_dist 0

        set valid_pairs 0

        # 计算12个距离

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

            set resB [expr {$i + 12}] ;# 链B的对应残基

            # 选择对应残基上的磷原子

            set selA [atomselect $mol "resid $i and element P" frame $fr]

            set selB [atomselect $mol "resid $resB and element P" frame $fr]

            # 检查选择是否有效

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

                # 安全获取坐标

                set coordsA [$selA get {x y z}]

                set coordsB [$selB get {x y z}]

                if {[llength $coordsA] == 3 && [llength $coordsB] == 3} {

                    # 计算距离

                    set dist [vecdist $coordsA $coordsB]

                    set dist [format "%.5f" $dist]

                    lappend distances $dist

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

                    incr valid_pairs

                    # 输出当前距离

                    puts -nonewline $fout [format "%12.5f" $dist]

                } else {

                    puts -nonewline $fout [format "%12s" "N/A"]

                    lappend distances "N/A"

                }

            } else {

                puts -nonewline $fout [format "%12s" "N/A"]

                lappend distances "N/A"

            }

            # 清理选择

            $selA delete

            $selB delete

        }

        # 计算平均直径

        if {$valid_pairs > 0} {

            set avg_R [expr {$sum_dist / double($valid_pairs)}]

            set avg_R [format "%.5f" $avg_R]

            puts $fout [format "%12.5f" $avg_R]

        } else {

            puts $fout [format "%12s" "N/A"]

        }

        # 显示进度

        if {$fr % 10 == 0} {

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

            flush stdout

        }

    }

    # 关闭文件

    close $fout

    puts "分析完成! 结果已保存至: $outfile"

    return 0

}

# 执行计算

if {[catch {calc_dna_radius} errmsg]} {

    puts stderr "发生错误: $errmsg"

    return -1

}