source: utils/developer/tcl/animatesphere.tcl

#!/usr/bin/tclsh
#
# This scripts displays the rolling sphere used in the tesselation procedure
# to compute the molecular surface via VMD's graphics interface.
# It is used during debugging and all its parameters are thrown out by an
# error message in CandidateForTesselation.cpp.


proc animate_sphere { molid atomid1 atomid2 atomid3 rad x1 y1 z1 x2 y2 z2 {step 25} } {
        global SphereRadius
        global vector1
        global vector2
        global CircleCenter
        global angle
        global radius
        global M_PI
        set SphereRadius $rad

        # avoid some 'animate dup' bug
        set molid [molinfo $molid get id]
        # check molid
        if {[molinfo $molid get numframes] != 1} {
                error "Can't find $molid or more than 1 frame."
        }
        # step must be integer and greater 2
        if {$step != int($step)} {
                error "step must be integer."
        }
        if {$step < 2} {
                error "step should be greater than 2."
        }
        # make step-1 copies of current frame
        for {set i 1} {$i <$step} {incr i} {
                animate dup frame 0 $molid
        }

        # select both atoms
        incr atomid1 -1
        incr atomid2 -1
        incr atomid3 -1
        set selid1 [atomselect $molid "index $atomid1" frame 0]
        set selid2 [atomselect $molid "index $atomid2" frame 0]
        set selid3 [atomselect $molid "index $atomid3" frame 0]

        # calculate origin of rotation
        set Atom1 [$selid1 get {x y z}]
        puts "Atom1 is $Atom1."
        set Atom2 [$selid2 get {x y z}]
        puts "Atom2 is $Atom2."
        set Atom3 [$selid3 get {x y z}]
        puts "Atom3 is $Atom3."
        set CircleCenter [vecadd [vecscale 0.5 [expr $Atom1]] [vecscale 0.5 [expr $Atom2]]]
        set CircleCenter [vecadd [vecscale 0.5 [expr $Atom1]] [vecscale 0.5 [expr $Atom2]]]
        puts "CircleCenter is at $CircleCenter."
        
        # calculate normal
        set Normaltmp [vecsub [expr $Atom1] [expr $Atom2]]
        set CirclePlaneNormal [vecnorm $Normaltmp]
        puts "CirclePlaneNormal is at $CirclePlaneNormal."

        # calculate radius
        set SphereCenter [vecsub [list $x1 $y1 $z1] $CircleCenter]
        puts "SphereCenter is at $SphereCenter."
        set OtherSphereCenter [vecsub [list $x2 $y2 $z2] $CircleCenter]
        puts "OtherSphereCenter is at $OtherSphereCenter."
        set radius [veclength $SphereCenter]
        puts "Radius of circle is $radius."
        set vector1 [vecnorm $SphereCenter]
        puts "vector1 is $vector1."

        # calculate other normal vector on plane, pointing out of triangle
        set vector2 [veccross $CirclePlaneNormal $vector1]
        set helper [vecsub $CircleCenter [expr $Atom3]]
        if {[vecdot $helper $vector2] < -1e-10} {
                puts "vector2 points inwards, rescaling."
                set vector2 [vecscale -1. $vector2]
        }
        puts "vector2 is $vector2."

        # get angle
        set angle [expr acos([vecdot $SphereCenter $OtherSphereCenter]/[veclength $SphereCenter]/[veclength $OtherSphereCenter])]

        # make angle in right direction (rotating outward in the sense of vector2
        set SKP [vecdot $OtherSphereCenter $vector2]
        puts "SKP with SearchDirection is $SKP."
        if {$SKP < -1e-10} {
                set angle [expr 2.*$M_PI - $angle]
        }
        puts "angle is $angle."

        global vmd_frame
        trace variable vmd_frame([molinfo top]) w update_current_sphere
        display update
        #enumerate_atoms 0
        return
}       

proc update_current_sphere {name index op} { 
        global SphereRadius
        global radius
        global angle
        global vector1
        global vector2
        global CircleCenter
        # draw sphere
        set CurrentAngle [expr ($angle/[expr [molinfo $index get numframes]-1])*[molinfo $index get frame]]
        set center [vecadd [vecscale $radius [vecadd [vecscale [expr cos($CurrentAngle)] $vector1] [vecscale [expr sin($CurrentAngle)] $vector2]]] $CircleCenter]
        #puts "center is now at $center with CurrentAngle is $CurrentAngle and radius $SphereRadius."
        draw material Transparent
        draw color silver
        if {[draw exists 2] != 0} {
                draw replace 2
        }
        draw sphere $center radius [expr $SphereRadius] resolution 30
        return
}

proc animate_sphere_off {} {
        global vmd_frame
        trace vdelete vmd_frame([molinfo top]) w update_current_sphere
        draw delete all
        animate delete beg 1 end 25 skip 0 [molinfo top]
        return
}

#enumerate_atoms 0
#animate_sphere 0 15 6 3. 4.00145 1.4968 1.28282 1.97354 -2.53854 -1.81601