Scripting Biomolecular Simulations with Tcl part of CS590v

1. Scripting Biomolecular Simulations with Tclpart of CS590v Michael McLennan Senior Research Scientist Rosen Center for Advanced Computing, Purdue University mmclennan@purdue.edu

2. Biomolecular Simulation developed by the Theoretical and Computational Biophysics Group in the Beckman Institute for Advanced Science and Technology at the University of Illinois at Urbana-Champaign Download from http://www.ks.uiuc.edu/Research/vmd/

3. What is scripting? • Type in commands • Write little programs • Like Unix “shell” scripts

4. button .b –text “Hello, World!” pack .b –pady 8 Tcl: language Tk: widget toolkit What is Tcl? • Developed by John Ousterhout at UC Berkeley • Released in 1989, currently maintained as Open Source • Millions of users worldwide • Used in commercial CAD tools and products (TiVo!) • More info: http://www.tcl.tk/

5. Application: Ubiquitin Ubiquitin 660 atoms Carbon backbone drawn in green

6. Application: Ubiquitin Ubiquitin 660 atoms Carbon backbone drawn in green (tube method)

7. Application: Ubiquitin Ubiquitin 660 atoms = starting point = after simulation

8. Application: Ubiquitin Ubiquitin 660 atoms Movement after simulation: = small change = medium change = large change

9. Coloring script for comparing molecules • A little script to color according to atomic movements: proc tutorialcoloring {} { set mol1 [lindex [molinfo list] 0] set mol2 [lindex [molinfo list] 1] set sel0 [atomselect $mol1 "alpha and protein"] ;# crystal set sel1 [atomselect $mol2 "alpha and protein"] ;# simulation set mylist {} foreach v0 [$sel0 get {x y z}] v1 [$sel1 get {x y z}] { set dx [expr [lindex $v0 0] - [lindex $v1 0]] set dy [expr [lindex $v0 1] - [lindex $v1 1]] set dz [expr [lindex $v0 2] - [lindex $v1 2]] set disp [expr ($dx*$dx + $dy*$dy + $dz*$dz)] lappend mylist $disp } $sel0 set beta $mylist } Tcl code!

10. Exit tclsh Add options to commands Output written on stdout by default Command to write out a string Output written without trailing newline Like Unix shell, but speaks Tcl tclsh + widgets VMD molecular visualization “Hello, World!” in Tcl • Start up a Tcl application: unix> tclsh unix> wish unix> vmd or or For example: unix> tclsh % puts “Hello, World!” Hello, World! % unix> tclsh % puts “Hello, World!” Hello, World! % puts –nonewline “hello” hello% unix> tclsh % puts “Hello, World!” Hello, World! % puts –nonewline “hello” hello% % % exit unix>

11. command arg arg … Tcl command syntax • First word is the command name • Remaining arguments depend on command syntax • Use double-quotes (“”) to wrap up text strings • Any line starting with hash (#) is a comment • Use semicolon (;) to separate multiple commands on same line • Use backslash (\) to extend a command onto multiple lines puts –nonewline “Hello, World!” # this is a comment puts “Hello”; puts “World!” puts –nonewline \ “Hello, World!”

12. Tcl Documentation http://www.tcl.tk/doc/

13. set x “World” puts “Hello, $x!” èHello, World! Uses the value of variable named x set pi 3.14159 puts “Value of pi is $pi” èValue of pi is 3.14159 set cmd puts $cmd “Hello, World!” èHello, World! set var x set $var “Universe” puts “Hello, $x!” èHello, Universe! Value can be a number Value can be a command name Value can be a variable name Tcl Variables • Variables hold numbers, strings, and other values:

14. Why didn’t we get a number here? èCircumference is 2*3.14159*10 expr 2+2 è4 Use the expr command to do math expr 2*$pi*$r è62.8318 The right way to get circumference Oops! Still not quite right! èCircumference is expr 2*3.14159*10 puts “Circumference is [expr 2*$pi*$r]” Execute command and substitute result in its place Doing Math in Tcl • Mathematical Expressions: set pi 3.14159 set r 10 puts “Circumference is 2*$pi*$r” puts “Circumference is expr 2*$pi*$r” èCircumference is 62.8318

15. c b a c = Ö a2+b2 In C you might say… In Tcl you say… Tcl Versus Other Languages • Pythagorean Theorem expressed in Tcl: set a 3 set b 4 set c [expr sqrt($a*$a + $b*$b)] • Use ()’s for functions/grouping • Use C-like math functions x = x + 1; set x [expr $x + 1] x += 2; incr x 2 set x [expr $x+2] Y = sin(2*x/(x-1)); set y [expr sin(2*$x/($x-1))]

16. puts “Hello, World!” puts {Hello, World!} Both handle multi-line strings The difference is important when you have substitutions: set x “World” puts “Hello, $x!” èHello, World! set x “World” puts {Hello, $x!} èHello, $x! Curly braces prevent substitutions! Tcl Quoting Rules • Tcl also supports {} quotes: puts “Hello, World!” puts {Hello, World!} Both keep text together

17. Problem #1 set x 10 puts x Problem #2 set set set puts $set Problem #3 set pi 3.14159 set area {expr $pi*$r*$r} puts “Area is: $area” Quiz èx èset èArea is: expr $pi*$r*$r

18. Command Scripts • Save a series of commands in a command script file: File: hello.tcl # This is a command script puts –nonewline “What’s your name?” set name [gets stdin] puts “Hello, $name!” Load the script interactively: Run as a program: unix> tclsh % source hello.tcl What’s your name?Fred Hello, Fred! % source hello.tcl What’s your name? unix> tclsh hello.tcl What’s your name?Fred Hello, Fred! unix>

19. Programming Statements Conditionals: if {$x > 0} { statements } if {$x > 0} { statements } elseif{$x < 0} { statements } else { statements } switch -regexp $x { a.*z {statements} [0-9]+ {statements} foo – bar {statements} } Looping: while {$x != 0} { statements } for {set x 0} {$x < 10} {incr x} { statements } foreach x {a b c d e} { statements } break continue Break out of loop Go back to top of loop

20. set x {a b c d e} llength $x è5 set x “a b c d e” llength $x è5 lappend x f lappend x g h puts $x èa b c d e f g h lindex $x 0 èa lindex $x 1 èb lindex $x end èh Space-separated list of values Nothing special about quotes Add to list stored in variable x Extract an element from a list Lists of Values • Variables can hold lists of values

21. proc lreverse {list} { set rlist “” set max [expr [llength $list]-1] for {set i $max} {$i >= 0} {incr i -1} { set v [lindex $list $i] lappend rlist $v } return $rlist } set x {a b c d e f g} set rx [lreverse $x] proc lreverse {list} { set rlist “” set max [expr [llength $list]-1] for {set i $max} {$i >= 0} {incr i -1} { set v [lindex $list $i] lappend rlist $v } return $rlist } set x {a b c d e f g} set rx [lreverse $x] proc lreverse {list} { set rlist “” set max [expr [llength $list]-1] for {set i $max} {$i >= 0} {incr i -1} { set v [lindex $list $i] lappend rlist $v } return $rlist } set x {a b c d e f g} set rx [lreverse $x] proc lreverse {list} { set rlist “” set max [expr [llength $list]-1] for {set i $max} {$i >= 0} {incr i -1} { set v [lindex $list $i] lappend rlist $v } return $rlist } set x {a b c d e f g} set rx [lreverse $x] Procedure declaration Empty list name argument list Substitutions (innermost to outermost) i = max; i >= 0; i-- Get element from list Add to end of return list Return value from procedure call Functions and Procedures • Reverse a list:

22. molecules top molecule VMD Commands • VMD is just like tclsh, but with extra Tcl commands: • vmd > molinfo list è0 1 vmd >atomselect top all èatomselect0 vmd >atomselect0 num è660 vmd >atomselect0 set radius 10 vmd >atomselect0 set radius 2 vmd >atomselect top alpha èatomselect1 vmd >atomselect1 set radius 5 vmd >atomselect top hydrophobic èatomselect2 vmd >atomselect2 set radius 10

23. VMD Commands • Use variables to store selections: vmd >set s [atomselect top all] vmd >$s num è660 vmd >$s set radius 2 vmd >set s2 [atomselect top alpha] vmd >$s2 num è76 vmd >$s2 set radius 5

24. Coloring script for comparing molecules • A little script to color according to atomic movements: proc tutorialcoloring {} { set mol1 [lindex [molinfo list] 0] set mol2 [lindex [molinfo list] 1] set sel0 [atomselect $mol1 "alpha and protein"] ;# crystal set sel1 [atomselect $mol2 "alpha and protein"] ;# simulation set mylist {} foreach v0 [$sel0 get {x y z}] v1 [$sel1 get {x y z}] { set dx [expr [lindex $v0 0] - [lindex $v1 0]] set dy [expr [lindex $v0 1] - [lindex $v1 1]] set dz [expr [lindex $v0 2] - [lindex $v1 2]] set disp [expr ($dx*$dx + $dy*$dy + $dz*$dz)] lappend mylist $disp } $sel0 set beta $mylist }