proc highlight_bases {molecule_id} { # one way to pick the atoms in the base set base {(not (backbone or name "[CO][1-5]['*]"))} # get the list of residues set residues [lsort -unique [[atomselect $molecule_id "all"] get residue]] # for each residue, get the locations of the base foreach residue $residues { # get the center of mass of the base set baseatoms [atomselect $molecule_id "residue $residue and $base"] set com [measure center $baseatoms weight mass] # pyrimidines and purines are drawn differently # get the atoms as if this is a pyrimidine (get N9 and N7) set c1 [atomselect $molecule_id "residue $residue and name N9"] set c2 [atomselect $molecule_id "residue $residue and name C8"] set box_length 4.0 # if there isn't an N9, then this is a purine, if {[$c1 num] == 0} { # use C2 and C6 set c1 [atomselect $molecule_id "residue $residue and name N1"] set c2 [atomselect $molecule_id "residue $residue and name C6"] set box_length 3.5 } # get the coordinates lassign [$c1 get {x y z}] p1 lassign [$c2 get {x y z}] p2 # find the coordinates of the corners of the box set v1 [vecsub $p2 $p1] set perp [vecnorm [veccross [vecsub $p2 $p1] [vecsub $com $p1]]] set v2 [vecscale $box_length [vecnorm [veccross $perp $v1]]] set perp [vecscale $perp 0.3] ########### draw a box for the base # use the right color lassign [$c1 get resname] resname draw color [colorinfo category Resname $resname] # draw the top of the slab draw triangle $p1 $p2 [vecadd $p2 $v2] draw triangle $p1 [vecadd $p2 $v2] [vecadd $p1 $v2] # draw the bottom draw triangle [vecadd $perp $p1] [vecadd $perp $p2] \ [vecadd $perp $p2 $v2] draw triangle [vecadd $perp $p1] [vecadd $perp $p2 $v2] \ [vecadd $perp $p1 $v2] # and the long sides draw triangle $p1 [vecadd $p1 $perp] [vecadd $p1 $perp $v2] draw triangle $p1 [vecadd $p1 $perp $v2] [vecadd $p1 $v2] draw triangle $p2 [vecadd $p2 $perp] [vecadd $p2 $perp $v2] draw triangle $p2 [vecadd $p2 $perp $v2] [vecadd $p2 $v2] } }