Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
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