################################################################################
# analysis-capsid-waterdensity.tcl
#
# * Calculate the density of water inside viral capsid cavity and that outside
#   of the capsid by using the volume data obtained from executing script
#   "analysis-capsidvolume.tcl"
#
#   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              "waterflux-cavity)"; \
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 box		0
set shell 	1

if {$box} {
  set minmax {{-145.222686768 -143.865936279 -145.225921631} {146.435546875 143.101989746 147.117950439}}
  foreach {min max} $minmax {} 
  foreach {xmin ymin zmin} $min {}
  foreach {xmax ymax zmax} $max {}

  set t 0.0
  set selc "water and same residue as \
	(x > [expr $xmin-$t] and x < [expr $xmax+$t] and y > [expr $ymin-$t] and y < [expr $ymax+$t] and \
	z > [expr $zmin-$t] and z < [expr $zmax+$t] )"; 
	# puts "atom selection of host particles wanted: $selc"
  set selcomb [atomselect $molID $selc]
  set watnum1 [$selcomb num]; # puts "number of water atoms created in the domain: $watnum1"
  unset upproc_var_$selcomb
  puts $rout "$b $watnum1"

} elseif {$shell} {

# 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 vout    [open ${prefix}-lvcapsid.dat r];
gets $vout lVcavity
gets $vout lVexterior
close $vout

set rout    [open ${prefix}-watercavity.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"

set hselc "water";
set hsel [atomselect ${molID} $hselc]
set nt [$hsel num]; puts "total number of water atoms = $nt"
$hsel delete

puts " frame Ravg Rmin Rmax Thickness Ncavity Nshell Noutside Ntotal Dcavity Doutside (kg/m3)"
puts $rout " frame Ravg Rmin Rmax Thickness Ncavity Nshell Noutside Ntotal Dcavity Doutside (kg/m3)"

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 Vcavity [expr 4.0/3.0*$PI*$min*$min*$min]
    set Vcavity [lindex $lVcavity $i]
    set Vexterior [lindex $lVexterior $i]
    set Voutside [expr [lindex $wsize 0]*[lindex $wsize 1]*[lindex $wsize 2] - 4.0/3.0*$PI*$avg*$avg*$avg - $Vexterior]
    # set Voutside [expr 0.5*[lindex $wsize 0]*[lindex $wsize 1]*[lindex $wsize 2] - 4.0/3.0*$PI*$avg*$avg*$avg - $Vexterior]

    set hselc "water and same residue as (within $b of protein)";
    # puts "$remark selection for host particles: $hselc"
    set hsel [atomselect ${molID} $hselc frame $i]
    set nshell [$hsel num];
    set lshell [$hsel get index]
    set mlist [$hsel get mass]
    set mt [total_mass $mlist]
    $hsel delete

    set hselc "(x*x+y*y+z*z) < $ravg2 and water and not same residue as (within $b of protein) ";
    # set hselc "water and same residue as ( (not index [concat $lshell]) 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 Dinside [expr $mt/$Vcavity*pow(10,27)/$Na]
    $hsel delete

    set hselc "(x*x+y*y+z*z) >= $ravg2 and water and not same residue as (within $b of protein)";
    # set hselc "water and same residue as ( (not index [concat $lshell]) 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 Doutside [expr $mt/$Voutside*pow(10,27)/$Na]
    $hsel delete

    set nt [expr $ncavity+$nshell+$noutside]

    ### output radius results
    puts [format "%5d  %12.5f  %12.5f  %12.5f  %12.5f  %10d  %10d %10d  %10d   %12.5f  %12.5f" $i $avg $min $max [expr $max-$min] $ncavity $nshell $noutside $nt $Dinside $Doutside]
    puts $rout [format "%5d  %12.5f  %12.5f  %12.5f  %12.5f  %10d  %10d  %10d  %10d  %12.5f  %12.5f" $i $avg $min $max [expr $max-$min] $ncavity $nshell $noutside $nt $Dinside $Doutside]
}; # for

close $rout
}; # if box

exit