# This script uses a 2d delaunay triangulation to draw the surface of
# a lipid membrane.  Here is an example usage:
#
# source lipid_surface.tcl
#
# # first locate the two distinct layers of lipids (if necessary)
# set groups [lipid_pids "DLPE LPPC"]
#
# # next calculate the normal direction (if you don't know it)
# set normal [lipid_normal $groups]
#
# # now draw the lipid surfaces, removing triangles with edges
# # longer than 20A.
# set sel1 [atomselect top "index [lindex $groups 0]"]
# set sel2 [atomselect top "index [lindex $groups 1]"]
#
# lipid_clear
# graphics $graphicsmol color yellow
# lipid_surface $sel1 $normal 20
# lipid_surface $sel2 [vecscale $normal -1] 20
#
# Edit this variable to point to your copy of Hull! 
set hullpath /home/pgrayson/bin/hull

set found 0
foreach m [molinfo list] {
    if { [regexp "lipid surface" [molinfo $m get name]] } {
	set found 1
    }
}

if { ! $found } {
    set top [molinfo top]

    mol load graphics "lipid surface"

    material add membrane
    material change ambient membrane 0.000000
    material change specular membrane 1.000000
    material change diffuse membrane 0.580000
    material change shininess membrane 0.540000
    material change opacity membrane 0.620000

    molinfo $top set top 1
    display resetview
}

foreach m [molinfo list] {
    if { [regexp "lipid surface" [molinfo $m get name]] } {
	set graphicsmol $m
    }
}


# lipid_pids returns the ids of the P atoms in the two lipid layers,
# in two separate lists.
proc expand_pids { lipid_segs lipid_pids } {
    set sel [atomselect top \
      "(within 15 of index $lipid_pids) and name P and segname $lipid_segs"]
    return [$sel list]
}

proc lipid_pids { lipid_segs } {
    set lipid_pids ""

    set sel [atomselect top "name P and segname $lipid_segs"]
    set lipid_pids [lindex [$sel list] 0]

    set size [llength $lipid_pids]
    set lipid_pids [expand_pids $lipid_segs $lipid_pids]
    while { [llength $lipid_pids] != $size } {
        set size [llength $lipid_pids]
        set lipid_pids [expand_pids $lipid_segs $lipid_pids]
    }
    
    set sel [atomselect top "segname $lipid_segs and name P and not index $lipid_pids"]
    set second_pids [$sel list]
    return "\"$lipid_pids\" \"$second_pids\""
}

# usage: lipid_normal $lipid_pids
# returns an approximate normal to the lipid membrane
proc lipid_normal { lipid_groups } {
    set one_layer [lindex $lipid_groups 0]
    set p_atoms [atomselect top "index $one_layer"]
    set all_lipid [atomselect top "same residue as index $one_layer"]

    set rp [measure center $p_atoms]
    set rl [measure center $all_lipid]
    return [vecnorm [vecsub $rp $rl]]
}

proc lipid_clear { } {
    global graphicsmol

    graphics $graphicsmol delete all
}

# Give me a surface of lipids, and I will plot it!
proc lipid_surface { lipidsel normal cutoff } {
    global graphicsmol
    global tcl_platform
    global hullpath

    set cutoff2 [expr $cutoff*$cutoff]

    graphics $graphicsmol material membrane

    # find the heads
    set heads $lipidsel
    set coords [ $heads get {x y z} ]
    set headsx [ $heads get x ]
    set headsy [ $heads get y ]
    set headsz [ $heads get z ]
    set num [llength [$heads list]]

    # generate a RH-coordinate system {perp1, perp2, normal}
    set dotx [vecdot $normal {1 0 0}]
    set doty [vecdot $normal {0 1 0}]
    if { abs($dotx) > abs($doty) } { # use x or y, whichever is more different
        set perp1 [veccross $normal {0 1 0}]
    } else {
        set perp1 [veccross $normal {1 0 0}]
    }
    set perp2 [veccross $normal $perp1]
    set perp1 [vecnorm $perp1]
    set perp2 [vecnorm $perp2]

    # project the points into 2d and zero out the normal array
    set points ""
    set i 0
    foreach r $coords  {
        set r1 [vecdot $perp1 $r]
        set r2 [vecdot $perp2 $r]
	append points "$r1 $r2\n"
	set pointnormal($i) {0 0 0}
	incr i
    }

    # compute the delaunay triangles
    set delaunay [exec $hullpath -m 20 -X /dev/null -d << $points]
    set delaunay [split $delaunay "\n"]
    
    # read out the triangles
    foreach line $delaunay {
  	if { ! [regexp "^\%" $line] } {
            if { [llength $line] == 3 } {
                # get the vertex coordinates
                set p0 [lindex $line 0]
                set p1 [lindex $line 1]
                set p2 [lindex $line 2]
                
                if { $p0 != -1 && $p1 != -1 && $p2 != -1 } {

                    set coords0 [lindex $coords $p0]
                    set coords1 [lindex $coords $p1]
                    set coords2 [lindex $coords $p2]

                    # get the edgelengths
                    set dist0 [veclength2 [vecsub $coords0 $coords1]]
                    set dist1 [veclength2 [vecsub $coords1 $coords2]]
                    set dist2 [veclength2 [vecsub $coords2 $coords0]]

                    if { $dist0 < $cutoff2 && $dist1 < $cutoff2 && $dist2 < $cutoff2 } {
                        
                        # add the normal
                        set norm [vecnorm [veccross [vecsub $coords0 $coords1] \
                                [vecsub $coords0 $coords2]]]
                        if { [vecdot $norm $normal] < 0 } {
                            set norm [vecinvert $norm]
                            set temp $p0
                            set p0 $p1
                            set p1 $temp
                            set temp $coords0
                            set coords0 $coords1
                            set coords1 $temp
                        }
                        set pointnormal($p0) [vecadd $pointnormal($p0) $norm]
                        set pointnormal($p1) [vecadd $pointnormal($p1) $norm]
                        set pointnormal($p2) [vecadd $pointnormal($p2) $norm]

                        lappend triangles "\"$coords0\" \"$coords1\" \"$coords2\""
                        lappend trianglepoints "$p0 $p1 $p2"
                    }
                }
            } else { puts "Bad line: $line" }
        }
    }

    # normalize the normals
    foreach point [array names pointnormal] {
	if { [veclength2 $pointnormal($point)] > 0 } {
	    set pointnormal($point) [vecnorm $pointnormal($point)]
	}
    }

    foreach triangle $triangles ids $trianglepoints {
	set norm0 $pointnormal([lindex $ids 0])
	set norm1 $pointnormal([lindex $ids 1])
	set norm2 $pointnormal([lindex $ids 2])
	graphics $graphicsmol trinorm "[lindex $triangle 0]" "[lindex $triangle 1]" "[lindex $triangle 2]" "$norm0" "$norm1" "$norm2"
    }
}