num-nt-contacts.tcl

source loadRigidBody.tcl

proc abs {a} {
    return $a > 0. ? $a : [expr { -1.*$a}]
}
proc wrapDiff {coords1 coords2} {
    set cell { 300. 300. 300. }
    set r 0.
    for {set i 0} { $i < 3} { incr i } {

        set x1 [lindex $coords1 $i]
        set x2 [lindex $coords2 $i]
        set x  [abs [expr { $x1 - $x2 }]]
        set l  [lindex $cell $i]
        while { $x > [expr { 0.5 *$l }]  } {
            set x [ expr {$x - $l} ]
        }
        set r [expr { $r + $x*$x}]
    }
    return $r
}

proc MeasureContactPBC {dna protein} {

    set count 0
    set cutoff 20
    set cutoff2 [expr { $cutoff * $cutoff } ]
    foreach coords1 $dna {
        foreach coords2 $protein {
            set r [wrapDiff $coords1 $coords2]
            if {$r < $cutoff2} {
                set count [ expr { $count + 1} ]
                break
            }
        }
    }
    return $count
} 

set sel1 [atomselect $IDs "name P0"]
set sel2 [atomselect $IDs "name P1"]
set pro  [atomselect $rbIDs "name CA"]
#set id3 [$pro molid]
#package require pbctools
#pbc set {200. 200. 200.} -all
#pbc wrap -center origin -molid 0 -all
#pbc wrap -centersel {0 0 0} -molid $id1 -all
#pbc wrap -centersel {0 0 0} -molid $id2 -all
set ch [open num-nt-contacts.dat w]

for {set f 0} {$f < [molinfo $IDs get numframes]} {incr f} {
    animate goto $f
    $sel1 frame $f
    $sel2 frame $f
    #set c [$pro get {x y z}]
    #set id1 [$sel1 molid]
    #puts $id1
    #set id2 [$sel2 molid]
    #puts $id2
    #set id3 [$pro molid]
    #pbc wrap -center origin -molid 0 -now
    #pbc wrap -center origin -molid 1 -now
    #pbc wrap -center origin -molid 1 -now
    lassign [measure contacts 10 $sel1 $pro] is js
    set c1 [llength [lsort -unique $is]]
    #set c1 [MeasureContact [$sel1 get {x y z}] [$pro get {x y z}]]
    lassign [measure contacts 10 $sel2 $pro] is js
    set c2 [llength [lsort -unique $is]]
    #set c2 [MeasureContact [$sel2 get {x y z}] [$pro get {x y z}]]
    puts $ch "$f $c1 $c2"
}
close $ch