################################################################################
# analysis-sasa.tcl
#
#   Calculate SASA for trajectory frames.
#
#   Yinglong Miao
#   Center for Cell and Virus Theory
#   Department of Chemistry
#   Indiana University, Bloomington, USA
#   yimiao@indiana.edu 
#
#   6/4/2008
#
################################################################################

proc total_mass {mlist} {
  # calculate the total mass of an atom selection
  set sum 0
  foreach mass $mlist {
    set sum [expr {$sum + $mass}]
  }
  return $sum
}

# 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
if {$nreplicas > 1} { 
  set dcdfile ${prefix}${i}.dcd 
} else {
  set dcdfile ${prefix}.dcd 
}

mol new ${molname}.psf
mol addfile ${molname}.pdb
# mol addfile     ${molname}-out.restart.1000000.coor type namdbin first 0 last -1 step 1 waitfor all
mol addfile     $dcdfile first 0 last -1 step 1 waitfor all
set nf [molinfo top get numframes]; puts "Total number of frames loaded: $nf"

set molID "top"

# set pselc "protein"
set pselc "protein backbone and noh"
puts "$remark Atom selection for protein: $pselc"

# superimpose the trajectory to a reference structure (frame 0) regarding psel
set psel0 [atomselect ${molID} $pselc frame 0]
set psel [atomselect ${molID} $pselc]

set np [$psel num]; puts "Number of selected protein atoms: $np"

set selall [atomselect ${molID} all]

set rout    [open ${prefix}-sasa.dat w];

set PI 3.14159
set Na [expr 6.022*pow(10,23)]; puts "Na = $Na"

set b	2.4; puts "protein boundary = $b"

puts " frame Ravg Rmin Rmax Thickness Interior Exterior Difference Total"
puts $rout " frame Ravg Rmin Rmax Thickness Interior Exterior Difference Total"

for {set i 0 } {$i < $nf } { incr i } {
    $psel frame $i
    $selall frame $i

    # fit protein backbone
    if {$i > 0} {
      set mat [measure fit $psel $psel0 weight mass]
      $selall move $mat
    }; # if i

    ### Determine the center of mass of the molecule and store the coordinates
    set com [measure center $psel weight mass]; # puts "com: $com"
    set avg 0
    set max 0
    set min 1000000
    set rmsr 0
    ### Determine the distance of the farthest atom from the center of mass
    # puts "Number of selected atoms: [$sel num]"
    foreach pos [$psel 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 avg [expr $avg/$np]

    # water box size
    set minmax [measure minmax $selall]; # puts "water box minmax: $minmax"
    foreach {wmin wmax} $minmax { break }
    set wsize [vecsub $wmax $wmin]

    set ravg2 [expr $avg*$avg];
    set rmin2 [expr $min*$min];
    set rmax2 [expr $max*$max];

    # set hselc "water and (not index [concat $lshell]) and same residue as ( (x*x+y*y+z*z) < $ravg2 and (not within $b of protein) )";
    set hselc "protein backbone and noh and (x*x+y*y+z*z) < $ravg2 ";
    # puts "$remark selection for host particles: $hselc"
    set hsel [atomselect ${molID} $hselc frame $i]
    set ncavity [$hsel num];
    set lcavity [$hsel get index]
    set mlist [$hsel get mass]
    set mt [total_mass $mlist]
    set sasa_in [measure sasa $b $hsel -points lpts]
    $hsel delete

    set hselc "protein backbone and noh and (not index [concat $lcavity]) and (x*x+y*y+z*z) > $ravg2 ";
    # puts "$remark selection for host particles: $hselc"
    set hsel [atomselect ${molID} $hselc frame $i]
    set noutside [$hsel num];
    set mlist [$hsel get mass]
    set mt [total_mass $mlist]
    set sasa_ex [measure sasa $b $hsel -points lpts]
    $hsel delete

    set sasat [expr $sasa_in+$sasa_ex]
    set sasad [expr $sasa_ex-$sasa_in]

    ### output radius results
    puts [format "%5d  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f	%12.5f	%12.5f" $i $avg $min $max [expr $max-$min] $sasa_in $sasa_ex $sasad $sasat]
    puts $rout [format "%5d  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f	%12.5f	%12.5f" $i $avg $min $max [expr $max-$min] $sasa_in $sasa_ex $sasad $sasat]
}; # for

close $rout

exit