# given the selection, get the CAs and draw a distance plot


#VMD  --- start of VMD description block
#Name:
#  CA distance plot
#Synopsis:
#  Given a selection, plot the CA-CA distance plot
#Version:
#  1.0
#Uses VMD version:
#  1.1
#Ease of use:
#  2
#Procedures:
# <li> vecdist v1 v2 -- computes the norm of v2 - v1 (which are vectors
#   of length 3)
# <li> ca_dist selection -- makes a new graphics molecule which is
#   the CA-CA distance plot of the given selection
#Description:
#   Given a selection, this script finds the CA atoms and constructs a
# matrix whose elements are the distances from the CA of residue(i) to
# that on residue(j).  A new graphics molecule is created with the name
# "CA_distance(molid)", where molid is the index of the molecule used
# by the selection.  The distance matrix is plotted by residue number
# in the order they appeared in the PDB (or PSF, ...) file, but only
# the selected residues are shown.  The colors are determined by the
# color scale and may be changed accordingly.
#Example output:
# <a href="br.dist.gif"><img src="br.dist.small.gif">image</a>
# of bacteriorhodopsin with the CA-CA  on the bottom
#Files: 
# <a href="ca-dist.vmd">ca-dist.vmd</a>
#Author: 
# Andrew Dalke &lt;dalke@ks.uiuc.edu&gt;
#\VMD  --- end of block

# Input: two vectors of length three (v1 and v2)
# Returns: ||v2-v1||
# This is almost twice the speed as "veclength [vecsub $v1 $v2]"
# and gives ~10% overall performance increase
proc vecdist {v1 v2} {
    lassign $v1 x1 x2 x3
    lassign $v2 y1 y2 y3
    return [expr sqrt(pow($x1-$y1,2) + pow($x2-$y2,2) + pow($x3-$y3, 2))]
}


# Input: a selection
# Does: finds the CAs in the selection then computes and draws the
#   CA-CA distance grid with colors based on the color scale
# Returns: the id of the new graphics molecule containing the grid

# In terms of timing, 1/3 is spent in the first big loop,
# and the 2/3 in the second loop is shared between the Tcl commands
# and the three 'graphics' commands.  This routine takes about
# 70 seconds (on a Crimson) for 234 residues (which is 27495
# distance calculations and 54756*2 triangles!).
proc ca_distance {main_sel} {
    # get the CA atoms from the selection
    set mol [$main_sel molindex]
    set sel [atomselect $mol "@$main_sel and name CA"]

    # find distances between each pair
    set coords [$sel get {x y z}]
    set max 0
    set list2 [$sel list]

    foreach atom1 $coords id1 [$sel list] {
	foreach atom2 $coords id2 $list2 {
	    set dist($id1,$id2) [vecdist $atom2 $atom1]
	    set dist($id2,$id1) $dist($id1,$id2)
	    set max [max $max $dist($id1,$id2)]
	}
	lvarpop list2
	lvarpop coords
    }


    # draw the pretty graphic
    puts "Distances calculated, now drawing the distance map ..."
    mol load graphics "CA_distance($mol)"
    set gmol [molinfo top]
    # turn material characteristics off
    graphics $gmol materials off
    # i1 and j1 are i+1 and j+1; this speeds up construction of x{01}{01}
    set i 0
    set i1 1
    # preset the scaling factor (31.95 instead of 32 ensures there will be
    # no color 66 (for the max color), which is transparent)
    set scale [expr 31.95 / ($max + 0.)]
    set list2 [$sel list]
    foreach id1 [$sel list] {
	set j 0
	set j1 1
	set x00 "$i $j 0"
	set x10 "$i1 $j 0"
	set x11 "$i1 $j1 0"
	set x01 "$i $j1 0"
	foreach id2 $list2 {
	    set col [expr int($scale * $dist($id1,$id2)) + 34]
	    graphics $gmol color $col
	    graphics $gmol triangle $x00 $x10 $x11
	    graphics $gmol triangle $x00 $x11 $x01
	    incr j
	    incr j1
	    set x00 $x01
	    set x10 $x11
	    set x11 "$i1 $j1 0"
	    set x01 "$i $j1 0"
	}
	incr i
	incr i1
    }
    # put some numbers down to give an idea of what is where
    set resids [$sel get resid]
    set num [llength $resids]
    set start [lindex $resids 0]
    set end [lindex $resids [expr $num - 1]]
    graphics $gmol color yellow
    graphics $gmol text "0 0 0" "$start,$start"
    graphics $gmol text "$num 0 0" "$end,$start"
    graphics $gmol text "0 $num 0" "$start,$end"
    graphics $gmol text "$num $num 0" "$end,$end"
    return $gmol
}