Command Input File

for Tutorial on Using BTCL and TCL Commands To Manipulate the Geometry of Two Helices

This example file contains links to the complete descriptions of the BTCL commands. The BTCL language.

The annotated output file.

#BIOSYM btcl 3
#
#  Input File For Discover Generated By Insight Version 2.3.5 
#  Date:	   Fri Apr 22 13:56:30 1994
#  Host Name: iris24  	Host Type: iris
#  (with minor editing for version 95.0)

# define procedure to
# fix the relative position of two helices so that their axes are parallel
# and the line joining their centroids is perpendicular to the axes
#
proc fixorientation {molecule1 ca1 molecule2 ca2 mindist axis1var \
                 point1var axis2var point2var} {
    upvar 1 $axis1var axis1
    upvar 1 $point1var point1
    upvar 1 $axis2var axis2
    upvar 1 $point2var point2

    # find axes for molecules 1 and 2 by finding the least square fit lines
    # through the backbone atoms

    molGeom get lsqLine axis1 $ca1
    molGeom get lsqLine axis2 $ca2

    # find the angle between two lines

    geometry ang12 angle $axis1 $axis2
    
    # find the cross product between two lines 
    # so the axis of rotation can be defined

    geometry vec1 vector $axis1
    geometry vec2 vector $axis2
    vector vec3 cross $vec1 $vec2

    # rotate molecule2 along the axis with origin at centroid of axis2
    # and direction of vector perpendicular to plane defined by axis1 and axis2
    # to make the two helices parallel

    geometry point1 point $axis1
    geometry point2 point $axis2
    geometry rotateaxis "line $point2 $vec3"
    
    vector ang12 negate $ang12

    # change between angles and degrees because
    # Molgeom command takes angles in degrees while geometry command 
    # takes angles in radians
    vector ang12 degree $ang12
    molGeom rotate $rotateaxis $ang12 $molecule2

    # find the new axis of molecule 2  
    molGeom get lsqLine axis2 $ca2
    geometry vec2 vector $axis2
    geometry point2 point $axis2

    # find the vector which passes centroid of axis1 and perpendicular to 2

    geometry vect vector $axis2 $point1
    geometry perline "line $point1 $vect"
    geometry dist distance $point2 $perline
    vector dist negate $dist

    # translate the molecule along its own axis so 
    # two centroids of the axes will be the shortest distance apart

    molGeom translate $vec2 $dist $molecule2

    #
    # find the distance between the two helices and put them at $mindist apart
    #

    geometry point2 point $axis2
    geometry direction vector $axis1 $point2
    geometry dist distance $point1 $axis2

    vector dist subtract $mindist $dist

    molGeom translate $direction $dist $molecule2

    # find the new axis of molecule 2
    molGeom get lsqLine axis2 $ca2
    geometry point2 point $axis2
}

#   define procedure to
#   write out the heading for the table file
#
proc writeheading { tblfile } {
    puts $tblfile "!BIOSYM output_table 1"
    puts $tblfile "#"
    puts $tblfile "TITLE: Frame Number"
    puts $tblfile "MEASUREMENT TYPE: Number"
    puts $tblfile "UNITS OF MEASUREMENT: none"
    puts $tblfile "FUNCTION: Frame_number"
    puts $tblfile "#"
    puts $tblfile "TITLE: Rotate angle"
    puts $tblfile "MEASUREMENT TYPE: Angle"
    puts $tblfile "UNITS OF MEASUREMENT: radian"
    puts $tblfile "FUNCTION: Rotate_angle"
    puts $tblfile "#"
    puts $tblfile "TITLE: Spin angle"
    puts $tblfile "MEASUREMENT TYPE: Angle"
    puts $tblfile "UNITS OF MEASUREMENT: radian"
    puts $tblfile "FUNCTION: Spin_angle"
    puts $tblfile "#"
    puts $tblfile "TITLE: Separation distance"
    puts $tblfile "MEASUREMENT TYPE: Separation"
    puts $tblfile "UNITS OF MEASUREMENT: Angstroms"
    puts $tblfile "FUNCTION: Separation"
    puts $tblfile "#"
    puts $tblfile "TITLE: Total energy"
    puts $tblfile "MEASUREMENT TYPE: Energy"
    puts $tblfile "UNITS OF MEASUREMENT: kcal/mol"
    puts $tblfile "FUNCTION: Total energy"
    puts $tblfile "#" 
    puts $tblfile "#" 
}

#
#   MAIN COMMANDS
#
begin
#
#   set the maximum and minimum distances between two helices
# 
$mindist = 10.0
$maxdist = 15.0
$nummove = 5
$moveinc = ( $maxdist - $mindist ) / $nummove
#
#   set the start, end, and angle increases for self spin of helix 2
#
$startspin = 0
$endspin = 360.0
$numspin = 4
$spininc = ( $endspin - $startspin ) / $numspin
#
#   set the start, end, and angle increases for rotation of helix 2 around
#   helix 1
#
$startrotate = 0
$endrotate = 360
$numrotate = 4
$rotateinc = ( $endrotate - $startrotate ) / $numrotate
#  
#   initialize the frame number for distance between 2 helices
#
$nframe = 0
$distance = $mindist
# 
#   open a table file to be read into Insight interface
#
set tblfile [ open $PROJECT.tbl w ]
writeheading $tblfile
#
#   define names for the molecules and subsets of backbone atoms   
#
set molecule1 "HELIX_1:*:Atom;*"
set molecule2 "HELIX_2:*:Atom;*"
#
#   define and get the subsets of backbone atoms of the two helices
#
subset define Subset "HELIX_1:*:backbone1" "CA,N,C"
subset define Subset "HELIX_2:*:backbone2" "CA,N,C"
subset get helix1ca "HELIX_1:*:backbone1"
subset get helix2ca "HELIX_2:*:backbone2"

#
#   fix the orientations so two helix axes are parallel and centroids of
#   the helix are lined up
#

fixorientation $molecule1 $helix1ca $molecule2 $helix2ca $mindist \
    axis1 point1 axis2 point2
#
# print out the axes to check the results
#
echo $axis1
echo axis1 = [geometry axis1]
echo axis2 = [geometry axis2]
geometry dist distance $axis1 $axis2
echo line distance = [geometry dist]
geometry dist distance $point1 $point2
echo point distance = [geometry dist]

#   
#   write out the coordinates after reorientation
#
writeFile coordinate filename = .cor

#
#   change relative orientations of two helices and 
#   calculate energy at each orientation
#

#
#   rotate molecule 2 around molecule 1
#
for {$i = 0} {$i < $numrotate} {incr i} {
    $rotateang = $startrotate + $i * $rotateinc
    if {$i > 0} {
	molGeom rotate $axis1 $rotateinc $molecule2
    }
    #   find axis of molecule2 and the direction of translation again
    molGeom get lsqLine axis2 $helix2ca
    geometry point2 point $axis2
    geometry transdirect vector $axis1 $point2
    #
    #   rotate molecule 2 around its own axis
    #
    for {$j = 0} {$j <= $numspin} {incr j} {
	$spinang = $startspin + $j * $spininc
	if {$j > 0} {
	    molGeom rotate $axis2 $spininc $molecule2
	}
	#   move molecule2 so it is at minimum distance from molecule 1 again
	$movedist =  $mindist - $distance
	molGeom translate $transdirect $movedist $molecule2
	#
	#
	#   move molecule 2 along the direction defined by the perpendicular 
	#   direction between axis 1 and axis 2
	#
	for {$k = 0} {$k <= $nummove} {incr k} {

	    if {$k > 0} {
		molGeom translate $transdirect $moveinc $molecule2
	    }
	    #   find the new axis of molecule2 again
	    molGeom get lsqLine axis2 $helix2ca
	    #   find the distances between two axes
	    geometry dist distance $axis1 $axis2
	    set distance [object dist]
            if {$j == $numspin} continue
	    #   evaluate energy
	    set energyval [energy]
	    #   write the coordinates in an archive file
	    $nframe = $nframe + 1
	    writeFile archive frame = $nframe filename = .arc
	    #   write the heading in output file 
	    if {$nframe == 1} {
		puts [ format " \n " ]
		puts [ format " %12s %12s %12s %12s %12s " \
			   Frame Rotate Spin Distance Energy ]
		puts [ format " \n " ]
	    }
	    #   write the energy into the output file
	    puts [ format " %12d %12.2f %12.2f %12.2f %12.6g " \
		       $nframe $rotateang $spinang $distance $energyval ]
	    #   write the energy into the table file for sorting and graphing
	    puts $tblfile "$nframe $rotateang $spinang \
                                    $distance $energyval"
	} 
    }
}

Main access page

Advanced-Use (BTCL) access

BTCL - Tutorial Access

Helix Lesson, Introduction Helix Lesson, Lesson