Center of Mass of Retinal

  VMD easily animates the retinal being pulled from the bR.  There is
more you might want to do with this information.  You may want to list
the displayment through time of the center of mass of the retinal from
its first position.  Here is how to do that.

proc measure_com_displacement {sel} {
	# save the current frame number
	set old_frame [$sel frame]

	# how many frames are available?
	set molid [$sel molid]
	set num_frames [molinfo $molid get numframes]

	# go to the first frame
	$sel frame 0
	# compute the center of mass
	set com0 [measure center $sel weight mass]

	# loop over each frame
	for {set i 0} {$i < $num_frames} {incr i} {
		# go to that frame
		$sel frame $i
		# and compute the center of mass
		set com [measure center $sel weight mass]

		# figure the distance
		set dist [veclength [vecsub $com $com0]]
		# and print it out
		puts "timestep $i displacement: $dist"
	}

	# restore the original frame number
	$sel frame $old_frame
}


# here's an example use with the bR DCD file
set retinal_fragment [atomselect top \
    {resid 216 and not same fragment as index 3369}]
measure_com_displacement $retinal_fragment


   This is the output

timestep 0 displacement: 0.0
timestep 1 displacement: 4.00543e-05
timestep 2 displacement: 0.298418
   [...]



				Displaying the path

Here's a related script.  This will show the path of the center of
mass of the selection over time as a line.  The code is almost the
same as the previous example except that it creates a graphics
molecule and it draws a line instead of computing the distance.  It
also takes an additional parameter to specify the color of the path.
By default it is red.

proc plot_com_path {sel {color red}} {
	# create a new graphics molecule
	mol load graphics "com path of [$sel text]"
	set graphics [molinfo top]
	# set the default color
	graphics $graphics color $color

	# save the current frame number
	set old_frame [$sel frame]

	# how many frames are available?
	set molid [$sel molid]
	set num_frames [molinfo $molid get numframes]

	# go to the first frame
	$sel frame 0
	# compute the center of mass
	set prev [measure center $sel weight mass]

	# loop over each frame, except the first
	for {set i 1} {$i < $num_frames} {incr i} {
		# go to that frame
		$sel frame $i
		# and compute the center of mass
		set current [measure center $sel weight mass]

		# draw a line from the previous position to the current one
		graphics $graphics line $prev $current
		# the current position becomes the previous one
		set prev $current
	}

	# restore the original frame number
	$sel frame $old_frame
}


# Again, an example of use for retinal is to watch the center of
# mass of the retinal fragment in red
plot_com_path $retinal_fragment

# Also, watch the path of the C3 and the C15
set C3  [atomselect [$retinal_fragment molid] "resid 216 and name C3"]
set C15 [atomselect [$retinal_fragment molid] "resid 216 and name C15"]
# C3 in yellow, C15 in green
plot_com_path $C3  yellow
plot_com_path $C15 green