################################################################################
# analysis-56mers-TR.tcl
#
#   Analyize the rigid-body translation and rotation of capsomeres in CCMV capsid
#
#   Yinglong Miao
#   Center for Cell and Virus Theory
#   Department of Chemistry
#   Indiana University, Bloomington, USA
#   yimiao@indiana.edu 
#
#   2/25/2009
#
################################################################################
# set OP calculation parameters
set remark              "analysis)"; \
puts "$remark --------------------------------------------------------------------------------"

set parafile            molOP.dat;
set fpara               [open $parafile r]; \
puts "$remark Open file $parafile for reading parameters"
gets $fpara             molname; puts "$remark molecule name prefix: $molname"
gets $fpara             line; puts "$remark suffix line: $line"
set suffix	  	[lindex $line 0]; puts "$remark molecule name suffix: $suffix"
set nreplicas     	[lindex $line 1]; puts "$remark Number of replica MD trajectories = $nreplicas."
set dcdstep     	[lindex $line 2]; puts "$remark Frame step for reading DCD trajectory: $dcdstep"
close $fpara

set prefix        ${molname}-${suffix}; puts "file prefix: $prefix"
set logfile             vmd.log; puts "$remark logfile: $logfile"

# load molecules
mol new         ${molname}.psf waitfor all 
mol addfile     ${molname}.pdb
# mol addfile     1cwp-capsid56-vmd-H.pdb
# mol addfile     ${molname}-out-step0.1-1ns.dcd first 1 last -1 step 2 waitfor all
mol addfile     ${prefix}.dcd first 0 last -1 step $dcdstep waitfor all
set nf [molinfo top get numframes]; puts "Total number of frames loaded: $nf"

################################################################################
# atom selection
# set selection "name CA"; puts "atom selection: $selection"
# set selection "protein and backbone and noh"; puts "atom selection: $selection"

# obtain segid list of pentamers and hexamers
set np	12
set nh	20
set nps	180
# create segid list for pentamers
set pentamers {}
for {set i 0 } {$i < $np } { incr i } {
  set p [list [expr $i*5+1] [expr $i*5+2] [expr $i*5+3] [expr $i*5+4] [expr $i*5+5]]
  lappend pentamers $p
}
# puts "pentamers: $pentamers"

# create segid list for hexamers
set hexamers {}
for {set i 0 } {$i < $nh } { incr i } {
  set h [list [expr $np*5+$i*6+1] [expr $np*5+$i*6+2] [expr $np*5+$i*6+3] [expr $np*5+$i*6+4] [expr $np*5+$i*6+5] [expr $np*5+$i*6+6]]
  lappend hexamers $h
}
# puts "hexamers: $hexamers"

# open output files
set rout    [open ${prefix}-transform.dat w];
set r5out    [open ${prefix}-5mersT.dat w];
set r6out    [open ${prefix}-6mersT.dat w];

set selection "protein and backbone and noh"; # puts "atom selection: $selection"
set selall [atomselect top $selection]; # puts ""
set allframe0 [atomselect top $selection frame 0]; # puts ""

set allcen0 [measure center $allframe0 weight mass]; # puts "system COM: $allcen0"

set Pi	[expr acos(-1)]; 
# frame loop
puts "frame	dl5avg	dang5avg	dl6avg	dang6avg"
puts $rout "frame	dl5avg	dang5avg	dl6avg	dang6avg"

set l5	{}
set l5t	{}
set l5r	{}

set l6	{}
set l6t	{}
set l6r	{}

for {set j 1 } {$j <= $np } { incr j } {
  lappend l5 P${j}
}
for {set j 1 } {$j <= $nh } { incr j } {
  lappend l6 H${j}
}
puts $r5out "frame	[concat $l5]"
puts $r6out "frame	[concat $l6]"

for {set i 1 } {$i < $nf } { incr i } {
    # puts "frame: $i"
    # ----------------------------------------------------------
    # fit the entire structure 
    $selall frame $i
    $selall move [measure fit $selall $allframe0]

    # ----------------------------------------------------------
    # compute averaged capsomere translation distances
    # pentamer
    set lenl {}
    set dangl {}
    # set danglall {}
    for {set j 0 } {$j < $np } { incr j } {
	set selection "protein and backbone and noh and segid [lindex $pentamers $j]"; # puts "atom selection: $selection"
	set frame0 [atomselect top $selection frame 0]; # puts ""
	set cen0 [measure center $frame0 weight mass]; # puts "selected component COM of frame 0: $cen0"

	set sel [atomselect top $selection]; # puts ""
        $sel frame $i
        set cen [measure center $sel weight mass]; # puts "selected component COM of frame $i: $cen"

    	# ----------------------------------------------------------
    	# compute capsomere (l,theta,phi)
    	set v [vecsub $cen $allcen0];
    	set l [veclength $v]; 
    	# set dotl [expr [vecdot $cen $cen0]/($l*$l0)];
    	foreach {x y z} $v { break }
    	if {$x >= 0} {
      	  set theta [expr atan($y/$x)]
    	} else {
      	  set theta [expr atan($y/$x)+$Pi]
    	}; # if x
    	set phi [expr acos($z/$l)];
    	# ----------------------------------------------------------
    	# compute the translation distance
    	set dl [veclength [vecsub $cen $cen0]]; 
    	# ----------------------------------------------------------
    	# compute the rotation angle
    	set R1 [transaxis z [expr -$theta] rad]
    	set R2 [transaxis y [expr -$phi] rad]
    	set R4 [transaxis y $phi rad]
    	set R5 [transaxis z $theta rad]
    	# ----------------------------------------------------------
    	# measure the transformation matrix for $sel
    	set mat [measure fit $sel $frame0]; 
    	set R3 [transmult $R2 $R1 $mat $R5 $R4]
    	set dangc [expr acos([lindex [lindex $R3 0] 0])]
    	set dangs [expr -asin([lindex [lindex $R3 1] 0])]
    	set dang [expr ($dangc+$dangs)/2*180/$Pi];
    	# puts "$i	$l	[expr $theta*180/$Pi]	[expr $phi*180/$Pi]	$dl	[expr $dang*180/$Pi]"
    	puts "pentamer: $j	$l	[expr $theta*180/$Pi]	[expr $phi*180/$Pi]	$dl	[expr $dangc*180/$Pi]	[expr $dangs*180/$Pi]"

	# save computed numbers
	lappend lenl $dl
	if { $dang != 0 } { lappend dangl $dang }

	$frame0 delete
	$sel delete
    }; # for 
    set ndl5	[llength $lenl]
    set ndang5	[llength $dangl]
    if {$ndl5} { set dl5avg [expr [vecsum $lenl]/$ndl5] }
    if {$ndang5} { set dang5avg [expr [vecsum $dangl]/$ndang5] }
    puts $r5out	"$i	$lenl"
    unset lenl dangl
    # puts "ndl5 = $ndl5	ndang5 = $ndang5"
    # hexamer
    set lenl {}
    set dangl {}
    for {set j 0 } {$j < $nh } { incr j } {
	set selection "protein and backbone and noh and segid [lindex $hexamers $j]"; # puts "atom selection: $selection"
	set frame0 [atomselect top $selection frame 0]; # puts ""
	set cen0 [measure center $frame0 weight mass]; # puts "selected component COM of frame 0: $cen0"

	set sel [atomselect top $selection]; # puts ""
        $sel frame $i
        set cen [measure center $sel weight mass]; # puts "selected component COM of frame $i: $cen"

    	# ----------------------------------------------------------
    	# compute capsomere (l,theta,phi)
    	set v [vecsub $cen $allcen0];
    	set l [veclength $v]; 
    	# set dotl [expr [vecdot $cen $cen0]/($l*$l0)];
    	foreach {x y z} $v { break }
    	if {$x >= 0} {
      	  set theta [expr atan($y/$x)]
    	} else {
      	  set theta [expr atan($y/$x)+$Pi]
    	}; # if x
    	set phi [expr acos($z/$l)];
    	# ----------------------------------------------------------
    	# compute the translation distance
    	set dl [veclength [vecsub $cen $cen0]]; 
    	# ----------------------------------------------------------
    	# compute the rotation angle
    	set R1 [transaxis z [expr -$theta] rad]
    	set R2 [transaxis y [expr -$phi] rad]
    	set R4 [transaxis y $phi rad]
    	set R5 [transaxis z $theta rad]
    	# ----------------------------------------------------------
    	# measure the transformation matrix for $sel
    	set mat [measure fit $sel $frame0]; 
    	set R3 [transmult $R2 $R1 $mat $R5 $R4]
    	set dangc [expr acos([lindex [lindex $R3 0] 0])]
    	set dangs [expr -asin([lindex [lindex $R3 1] 0])]
    	set dang [expr ($dangc+$dangs)/2*180/$Pi];
    	# puts "$i	$l	[expr $theta*180/$Pi]	[expr $phi*180/$Pi]	$dl	[expr $dang*180/$Pi]"
    	puts "hexamer: $j	$l	[expr $theta*180/$Pi]	[expr $phi*180/$Pi]	$dl	[expr $dangc*180/$Pi]	[expr $dangs*180/$Pi]"

	# save computed numbers
	lappend lenl $dl
	if {$dang != 0} { lappend dangl $dang }

	$frame0 delete
	$sel delete
    }; # for 
    set ndl6	[llength $lenl]
    set ndang6	[llength $dangl]
    if {$ndl6} { set dl6avg [expr [vecsum $lenl]/$ndl6] }
    if {$ndang6} { set dang6avg [expr [vecsum $dangl]/$ndang6] }
    puts $r6out	"$i	$lenl"
    unset lenl dangl
    # puts "ndl6 = $ndl6	ndang6 = $ndang6"

    puts "$i	$dl5avg	$dang5avg	$dl6avg	$dang6avg"
    puts $rout "$i	$dl5avg	$dang5avg	$dl6avg	$dang6avg"

}; # for i nf

# delete the molecule
mol delete top

# close output files
close $rout
close $r5out
close $r6out

exit

# old, not smart, algorighm
# computate averaged translation distances of capsomeres via averaging those of protein subunits

if {0} {

set capsid	0
set pentamer	0
set hexamer	0
set protomer	0
set hexamers	0
set protomers	0

set segs {}

if {$capsid} {
set segstart	1
set segend	180
set segstep	1
for {set i $segstart } {$i <= $segend } { incr i $segstep} {
  lappend segs $i
}; # for i
set prefix	${prefix}-cap
}; # if capsid

if {$pentamer} {
set segstart	31
set segend	35
set segstep	1
for {set i $segstart } {$i <= $segend } { incr i $segstep} {
  lappend segs $i
}; # for i
set prefix	${prefix}-5mer
}; # if pentamer

if {$hexamer} {
set segstart	151
set segend	156
set segstep	1
for {set i $segstart } {$i <= $segend } { incr i $segstep} {
  lappend segs $i
}; # for i
set prefix	${prefix}-6mer
}; # if hexamer

if {$protomer} {
  set segs {51 159 156}
  set prefix	${prefix}-pro
}; # if pentamer

if {$pentamers} {
set segstart	1
set segend	60
set segstep	1
for {set i $segstart } {$i <= $segend } { incr i $segstep} {
  lappend segs $i
}; # for i
set prefix	${prefix}-5mers
}; # if pentamers

if {$hexamers} {
set segstart	61
set segend	180
set segstep	1
for {set i $segstart } {$i <= $segend } { incr i $segstep} {
  lappend segs $i
}; # for i
set prefix	${prefix}-6mers
}; # if hexamers

set selection "protein and backbone and noh and segid [concat $segs]"; puts "atom selection: $selection"
set sel [atomselect top $selection]; puts ""
set frame0 [atomselect top $selection frame 0]; puts ""
set num [$sel num]; puts "Number of selected atoms: $num"

# frame loop
for {set i 1 } {$i <= $nf } { incr i } {
    # puts "frame: $i"
    # ----------------------------------------------------------
    # fit the entire structure 
    $selall frame [expr $i-1]
    $selall move [measure fit $selall $allframe0]
    # ----------------------------------------------------------
    # compute the translation distance
    # pentamer
    set lenl {}
    for {set seg 1} {$seg <= 60} {incr seg} {
	set selection "protein and backbone and noh and segid $seg"; # puts "atom selection: $selection"
	set frame0 [atomselect top $selection frame 0]; # puts ""
	set cen0 [measure center $frame0 weight mass]; # puts "selected component COM of frame 0: $cen0"

	set sel [atomselect top $selection]; # puts ""
        $sel frame [expr $i-1]
        set cen [measure center $sel weight mass]; # puts "selected component COM of frame $i: $cen"

	lappend lenl [veclength [vecsub $cen $cen0]]

	$frame0 delete
	$sel delete
    }; # for 
    set dl5avg [expr [vecsum $lenl]/60]
    unset lenl
    # hexamer
    set lenl {}
    for {set seg 61} {$seg <= 180} {incr seg} {
	set selection "protein and backbone and noh and segid $seg"; # puts "atom selection: $selection"
	set frame0 [atomselect top $selection frame 0]; # puts ""
	set cen0 [measure center $frame0 weight mass]; # puts "selected component COM of frame 0: $cen0"

	set sel [atomselect top $selection]; # puts ""
        $sel frame [expr $i-1]
        set cen [measure center $sel weight mass]; # puts "selected component COM of frame $i: $cen"

	lappend lenl [veclength [vecsub $cen $cen0]]

	$frame0 delete
	$sel delete
    }; # for 
    set dl6avg [expr [vecsum $lenl]/120]
    unset lenl
    puts "$i	$dl5avg	$dl6avg"
}; # for i nf

}; # if capsomeres