################################################################################
# analysis-capsid-cRadius-rmsd.tcl
#
# * Analyize the changes in the radii of viral capsid backbone, i.e., the
#   average, maximum and minimum distances between atoms in the capsid backbone
#   and the capsid center of mass, in a simulation trajectory
#
# * Calculate the RMSD between viral capsid snapthots along simulation
#   trajectory and a reference structure
#
#   Yinglong Miao
#   Center for Cell and Virus Theory
#   Department of Chemistry
#   Indiana University, Bloomington, USA
#   yimiao@indiana.edu 
#
#   5/20/2008
#
################################################################################
# set OP calculation parameters
set remark              "analysis-Radius)"; \
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-NT0-vmd-H-wb-0.2M-2ns.pdb type pdb first 0 last -1 step 1 filebonds 1 autobonds 1 waitfor all
if {0} {
mol addfile 1cwp-swln1-capsid56-vmd-H.pdb type pdb first 0 last -1 step 1 filebonds 1 autobonds 1 waitfor all
mol addfile 1cwp-swln2-capsid56-vmd-H.pdb type pdb first 0 last -1 step 1 filebonds 1 autobonds 1 waitfor all
mol addfile 1cwp-capsid56-vmd-H.pdb type pdb first 0 last -1 step 1 filebonds 1 autobonds 1 waitfor all
}
# mol addfile     1cwp-capsid56-s5-vmd-H-out-MDOPX-step1-10ns.dcd first 0 last 8 step $dcdstep waitfor all
# 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"
# set selection "protein"; 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"

# open output files
set rout    [open ${prefix}-dRadius.dat w];
set rmsdout [open ${prefix}-rmsd.dat w];

puts $rout      "  t      Ravg(A)      Rmin(A)       Rmax(A)	cRavg(A)      cRmin(A)       cRmax(A)"
puts $rmsdout   "  t      RMSD(A)"

# frame loop
for {set i 1 } {$i <= $nf } { incr i } {
    puts "frame: $i"
    $sel frame [expr $i-1]
    # ----------------------------------------------------------
    # compute RMSD
    $sel move [measure fit $sel $frame0]
    puts $rmsdout [format "%5d  %12.5f" $i   [measure rmsd $sel $frame0]]
    # ----------------------------------------------------------
    # compute radius values
    ### Determine the center of mass of the molecule and store the coordinates
    set com [measure center $sel weight mass]
    set avg 0
    set max 0
    set min 1000000
    ### Determine the distance of the farthest atom from the center of mass
    # puts "Number of selected atoms: [$sel num]"
    foreach pos [$sel get {x y z}] {
      set dist [veclength [vecsub $pos $com]]
      set avg [expr $avg+$dist]
      if {$dist > $max} {set max $dist}
      if {$dist < $min} {set min $dist}
      set l2 [veclength2 $pos]
      }
    set avg [expr $avg/$num]
    ### output radius results
    if { $i == 1 } {
      set avg0	$avg
      set min0  $min
      set max0	$max
    }
    set davg	[expr $avg-$avg0]
    set dmin	[expr $min-$min0]
    set dmax	[expr $max-$max0]
    puts $rout [format "%5d  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f" $i   $avg    $min    $max	$davg	$dmin	$dmax]
}; # for i nf

# delete the molecule
mol delete top

# close output files
close $rout
close $rmsdout

exit