#!/usr/bin/perl

@temp=qw(280.00 286.04 292.18 298.41 304.74 311.17 317.69 324.31 331.03 337.85 344.79 351.82 358.97 366.22 373.59 381.07 388.67 396.40);

foreach $i (0..$#temp) {

  open OUT, ">pt-criteria.in";

  print OUT "parm ../../build/full.topo.hmr [traj]\n\nstrip :WAT,Na+,K+,Cl-\n\nautoimage origin\n\n";
  print OUT "rmsd :1,2&!\@H= first\n\n";

  print OUT "nativecontacts :1\@N1,C6,C5,C4,N4,N3,C2,O2 :2\@N1,C6,C5,C4,N4,N3,C2,O2 mindist\n";
  print OUT "distance COM :1\@N1,C6,H6,C5,H5,C4,N4,H41,H42,N3,C2,O2 :2\@N1,C6,H6,C5,H5,C4,N4,H41,H42,N3,C2,O2\n\n";

  print OUT "vector v_base2 :2\@N1,C6,C5,C4,N4,N3,C2,O2 corrplane\n";
  print OUT "vector v_base1 :1\@N1,C6,C5,C4,N4,N3,C2,O2 corrplane\n";
  print OUT "vectormath name normalangle vec1 v_base1 vec2 v_base2 dotangle\n\n";
  foreach $p (1..$ARGV[0]) {
    print OUT "trajin ../../traj.$p.01 remdtraj remdtrajtemp $temp[$i] trajnames ";
    foreach $k (1..$#temp) {
      $j=sprintf "%02.0f", $k+1;
      print OUT "../../traj.$p.$j,";
    }
   print OUT "\n\n";
  }

  print OUT "\nrun\n\n";

  print OUT "writedata criteria-raw_$temp[$i].dat Contacts_00001[mindist] COM normalangle\n";

  system "cpptraj -i pt-criteria.in > pt-criteria.log";

  open INP, "< criteria-raw_$temp[$i].dat";
  open OUT1, ">> temperature_dG_.dat";
  open OUT2, "> stacking_$temp[$i].dat";
  open OUT3, ">> temperature_stacked-Percentage.dat";

  @tot = 0;
  @stacked = 0;

  while (<INP>) {
          @col = split " ",$_;
          if ($col[0] ne '#Frame') {
                  $tot ++;
                  if (($col[1] < 4) && ($col[2] < 5) && (($col[3] < 45) || ($col[3] > 135)))  {
                          $stacked ++;
                  }
          }
  }
  $nonstacked = $tot - $stacked;
  $percentage = ($stacked / $tot) * 100;
  $relative = ($stacked/$nonstacked);
  if ($relative > 0) {
        $dG = -1.9872041 * 0.298 * log ($relative);
  }
  else {$dG = "inf";}
  print OUT2 "Total Frames = $tot\n";
  print OUT2 "Total Stacked = $stacked\n";
  print OUT2 "Total Non-Stacked = $nonstacked\n";
  print OUT2 "Stacked percetage = $percentage%\n";
  print OUT2 "Free energy = $dG \n";
  printf OUT1 "%3.2f    %3.2f\n", $temp[$i],$dG;
  printf OUT3 "%3.2f    %3.2f\n", $temp[$i],$percentage;
  close INP;
  close OUT2;
}
system "rm criteria*.dat";