################################################################################
# analysis-capsidvolume.tcl
#
# * Calculate the volume of a viral capsid and that inside the capsid cavity.
#   In case the capsid is solvated in a water box, the volume of the host medium
#   outside of the capsid is also calculated.
#
#   Yinglong Miao
#   Center for Cell and Virus Theory
#   Department of Chemistry
#   Indiana University, Bloomington, USA
#   yimiao@indiana.edu 
#
#   5/20/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              "capsid-volume)"; \
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"

# load combined structure
mol new /home/long/namd/vmd-scripts/combine.psf waitfor all
mol addfile /home/long/namd/vmd-scripts/combine.pdb waitfor all

# Set molID for loaded structures
set molIDs [molinfo list]; # puts "loaded molecule ID's: $molIDs"
set molID [lindex $molIDs 0]
set molIDcomb [lindex $molIDs 1]

# set molID "top"

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

# set selections for the capsid trajectory
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 PI 3.14159
set Na [expr 6.022*pow(10,23)]; puts "Na = $Na"

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

set pselc "protein"
set pselall [atomselect ${molID} $pselc]

# set selections for the combined structure
set pselcomb [atomselect ${molIDcomb} $pselc]

set hselc "water";
set hsel [atomselect ${molIDcomb} $hselc]
set nt [$hsel num]; puts "total number of water atoms in combine.psf = $nt"
set mlist [$hsel get mass]; puts ""
set mcomb [total_mass $mlist]
set Vcomb [expr 338.191009521*338.175003052*338.201004028]; puts "Vcomb = $Vcomb A^3"
set Dwater1 [expr $mcomb/$Vcomb*pow(10,27)/$Na] ; puts "water density = $Dwater1 kg/m3"
set Dwater [expr $mcomb/$Vcomb] ; puts "water density = $Dwater g/mol/A3"
$hsel delete

set rout    [open ${prefix}-capsidvolume.dat w];
set vout    [open ${prefix}-lvcapsid.dat w];

puts " frame Ravg Rmin Rmax Thickness Vinterior Vexterior Vcavity Vshell Vtotal (A^3) Nshell Ntotal"
puts $rout " frame Ravg Rmin Rmax Thickness Vinterior Vexterior Vcavity Vshell Vtotal (A^3) Nshell Ntotal"

set lVcavity {}
set lVexterior {}

for {set i 0 } {$i < $nf } { incr i } {
    $psel frame $i
    $pselall 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];

    $pselcomb set {x y z} [$pselall get {x y z}]

    set hselc "water and same residue as (within $b of protein)";
    # puts "$remark selection for host particles: $hselc"
    set hsel [atomselect ${molIDcomb} $hselc frame $i]
    set nshell [$hsel num];
    set lshell [$hsel get index]
    set mlist [$hsel get mass]; 
    set mt [total_mass $mlist]
    # set Dshell [expr $mt/$Vshell*pow(10,27)/$Na]
    set Vshell [expr $mt/$Dwater]
    $hsel delete

    # 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 "water and same residue as ((index [concat $lshell]) and (x*x+y*y+z*z) < $ravg2 )";
    # set hselc "water and same residue as (index [concat $lshell] and (x*x+y*y+z*z) < $ravg2) ";
    # set hselc "water and index [concat $lshell] and (x*x+y*y+z*z) < $ravg2 ";
    # set hselc "water and index [concat $lshell] and same residue as (within $b of protein and (x*x+y*y+z*z) < $ravg2 )";
    # set hselc "water and same residue as (within $b of protein and (x*x+y*y+z*z) < $ravg2 )";
    # set hselc "index [concat $lshell] and (x*x+y*y+z*z) < $ravg2";
    set hselc "(x*x+y*y+z*z) < $ravg2 and water and same residue as (within $b of protein)";
    # puts "$remark selection for host particles: $hselc"
    set hsel [atomselect ${molIDcomb} $hselc frame $i]
    set ninterior [$hsel num];
    set linterior [$hsel get index]
    set mlist [$hsel get mass]; 
    set mt [total_mass $mlist]
    set Vinterior [expr $mt/$Dwater]
    $hsel delete

    # set hselc "water and (not index [concat $lshell]) and same residue as ( (not index [concat $lshell]) and (x*x+y*y+z*z) > $ravg2 )";
    # set hselc "water and same residue as (index [concat $lshell] and (x*x+y*y+z*z) > $ravg2) ";
    # set hselc "water and index [concat $lshell] and (x*x+y*y+z*z) > $ravg2 ";
    # set hselc "water and index [concat $lshell] and same residue as (within $b of protein and (x*x+y*y+z*z) > $ravg2 )";
    # set hselc "index [concat $lshell] and (x*x+y*y+z*z) > $ravg2 ";
    set hselc "(x*x+y*y+z*z) > $ravg2 and water and same residue as (within $b of protein)";
    # puts "$remark selection for host particles: $hselc"
    set hsel [atomselect ${molIDcomb} $hselc frame $i]
    set nexterior [$hsel num];
    set mlist [$hsel get mass]; 
    set mt [total_mass $mlist]
    set Vexterior [expr $mt/$Dwater]
    $hsel delete

    # set Vcavity [expr 4.0/3.0*$PI*$min*$min*$min]
    # set Vshell [expr 4.0/3.0*$PI*($max*$max*$max-$min*$min*$min)]
    # set Vexterior [expr [lindex $wsize 0]*[lindex $wsize 1]*[lindex $wsize 2] - 4.0/3.0*$PI*$max*$max*$max]

    set Vcavity [expr 4.0/3.0*$PI*$avg*$avg*$avg-$Vinterior]
    set Vtotal [expr $Vinterior+$Vexterior]

    set ntotal [expr $ninterior+$nexterior]

    lappend lVcavity $Vcavity
    lappend lVexterior $Vexterior

    # puts [format "%5d  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f" $i $avg $min $max [expr $max-$min] $Vinterior $Vexterior $Vcavity $Vshell $Vtotal ]
    puts [format "%5d  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %10d  %10d" $i $avg $min $max [expr $max-$min] $Vinterior $Vexterior $Vcavity $Vshell $Vtotal $nshell $ntotal]
    puts $rout [format "%5d  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f  %12.5f" $i $avg $min $max [expr $max-$min] $Vinterior $Vexterior $Vcavity $Vshell $Vtotal $nshell $ntotal]

}; # for

puts $vout $lVcavity
puts $vout $lVexterior

close $rout
close $vout

exit