#!/bin/tcsh

# first argument is the number of runs

if ($1 == 1) then 
  rm -r analysis/
endif

mkdir analysis
mkdir analysis/rep-rmsd
mkdir analysis/remlog
mkdir analysis/conf-temp
mkdir analysis/stack-temp
mkdir analysis/suites-298
mkdir analysis/hbond-298
mkdir analysis/clusters-6

cd analysis

# RMSD per replica overlaps for convergence proof
cd rep-rmsd
cat << EOF > rep-rmsd.perl
#!/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){
  \$j=sprintf "%02.0f", \$i+1;

  open OUT, ">pt-rmsd.in";
  print OUT "parm ../../build/full.topo.hmr [traj]\nparm ../../../../../A-form.pdb [striped]\nreference ../../../../../A-form.pdb parm [striped] [a-form]\n\nstrip :WAT,Na+,K+,Cl-\n\nautoimage origin\n\nrmsd rms_2_aform :1,2&!\@H= out rmsd_\$j.dat ref [a-form]\n";
  close OUT;
  open OUT, ">>pt-rmsd.in";
  foreach \$p (1..\$ARGV[0]) {
    print OUT "trajin ../../traj.\$p.\$j \n";
  }
  system "cpptraj -i pt-rmsd.in > pt-rmsd.log";
}

system "rm pt-hist.in";
open HIST, ">>pt-hist.in";
foreach \$i (0..\$#temp){
  \$j=sprintf "%02.0f", \$i+1;
  print HIST "readdata rmsd_\$j.dat name replica_\$j\n";
}
print HIST "\nrun\n\n";
foreach \$i (0..\$#temp){
  \$j=sprintf "%02.0f", \$i+1;
  print HIST "hist replica_\$j,0,6,0.1 out hist_replicaRMSD_run$1.agr norm name replica_\$j\n";
}

system "cpptraj -i pt-hist.in > pt-hist.log";
EOF


#perl rep-rmsd.perl $1 &

rm rmsd_*.dat
cd ../

# remlog analysis
cd remlog

cat << EOF > remlog-analysis.perl
#!/usr/bin/perl

system "rm pt-remlog.in";

open OUT, ">>pt-remlog.in";
foreach \$p (1..\$ARGV[0]) {
    print OUT "readdata ../../remlog.\$p name REMLOG as remlog\n";
  }
print OUT "remlog REMLOG stats statsout roundtrips.dat reptime time_per_replica.dat out remlog_analysis.out repidx\n";

system "cpptraj -i pt-remlog.in > pt-remlog.log";
EOF

#perl remlog-analysis.perl $1
cd ../

# Conformational analysis in different temperatures
cd conf-temp

cat << EOF > ensemble.perl
#!/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);

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

print OUT "parm ../../build/full.topo.hmr [traj]\nparm ../../../../../A-form.pdb [striped]\nreference ../../../../../A-form.pdb parm [striped] [a-form]\n\nstrip :WAT,Na+,K+,Cl-\n\nautoimage origin\n\n";

print OUT "rmsd rms_2_aform :1,2&!\@H= ref [a-form] \n";
print OUT "distance dist_O2'-O4' :1\@O2' :2\@O4' \n";
print OUT "distance COM :1\@N1,C6,H6,C5,H5,C4,N4,H41,H42,N3,C2,O2 :2\@N1,C2,H2,N3,C4,C5,C6,N6,H61,H62,N7,C8,H8,N9 \n";

print OUT "dihedral g1 :1\@O5' :1\@C5' :1\@C4' :1\@C3'  type gamma  \n";
print OUT "dihedral d1 :1\@C5' :1\@C4' :1\@C3' :1\@O3'  type delta  \n";
print OUT "dihedral e1 :1\@C4' :1\@C3' :1\@O3' :2\@P    type epsilon\n";
print OUT "dihedral z1 :1\@C3' :1\@O3' :2\@P   :2\@O5'  type zeta   \n";
print OUT "dihedral c1 :1\@O4' :1\@C1' :1\@N1  :1\@C2   type chi    \n";
print OUT "dihedral a2 :1\@O3' :2\@P   :2\@O5' :2\@C5'  type alpha  \n";
print OUT "dihedral b2 :2\@P   :2\@O5' :2\@C5' :2\@C4'  type beta   \n";
print OUT "dihedral g2 :2\@O5' :2\@C5' :2\@C4' :2\@C3'  type gamma  \n";
print OUT "dihedral d2 :2\@C5' :2\@C4' :2\@C3' :2\@O3'  type delta  \n";
print OUT "dihedral c2 :2\@O4' :2\@C1' :2\@N9  :2\@C4   type chi    \n";

foreach \$nt (0..\$#temp) {
  system "mkdir \$temp[\$nt]";
}

foreach \$p (1..\$ARGV[0]) {
  print OUT "ensemble ../../traj.\$p.01 trajnames ";
  foreach \$i (1..\$#temp) {
    \$j=sprintf "%02.0f", \$i+1;
    print OUT "../../traj.\$p.\$j,";
  }
  print OUT "\n\n";
}

print OUT "\nrun\n\n";

foreach \$i (0..\$#temp) {
  print OUT "hist rms_2_aform%\$i,0,6,0.1 out rmsd-hist_Run\$ARGV[0].agr name rmsd_\$temp[\$i]K norm\n";
  print OUT "hist dist_O2'-O4'%\$i,0,10,0.1 out distO2O4-hist_Run\$ARGV[0].agr name Oxygen_distance_\$temp[\$i]K norm\n";
  print OUT "hist COM%\$i,0,10,0.1 out COM-hist_Run\$ARGV[0].agr name COM_distance_\$temp[\$i]K norm\n";

  print OUT "hist g1%\$i,-180,180,1 out dihed_g1-hist_Run\$ARGV[0].agr name gamma1_\$temp[\$i]K norm\n";
  print OUT "hist d1%\$i,-180,180,1 out dihed_d1-hist_Run\$ARGV[0].agr name delta1_\$temp[\$i]K norm\n";
  print OUT "hist e1%\$i,-180,180,1 out dihed_e1-hist_Run\$ARGV[0].agr name epsilon1_\$temp[\$i]K norm\n";
  print OUT "hist z1%\$i,-180,180,1 out dihed_z1-hist_Run\$ARGV[0].agr name zeta1_\$temp[\$i]K norm\n";
  print OUT "hist c1%\$i,-180,180,1 out dihed_c1-hist_Run\$ARGV[0].agr name chi1_\$temp[\$i]K norm\n";
  print OUT "hist a2%\$i,-180,180,1 out dihed_a2-hist_Run\$ARGV[0].agr name alpha2_\$temp[\$i]K norm\n";
  print OUT "hist b2%\$i,-180,180,1 out dihed_b2-hist_Run\$ARGV[0].agr name beta2_\$temp[\$i]K norm\n";
  print OUT "hist g2%\$i,-180,180,1 out dihed_g2-hist_Run\$ARGV[0].agr name gamma2_\$temp[\$i]K norm\n";
  print OUT "hist d2%\$i,-180,180,1 out dihed_d2-hist_Run\$ARGV[0].agr name delta2_\$temp[\$i]K norm\n";
  print OUT "hist c2%\$i,-180,180,1 out dihed_c2-hist_Run\$ARGV[0].agr name chi2_\$temp[\$i]K norm\n\n";
 
  print OUT "hist rms_2_aform%\$i,0,6,0.1 out ./\$temp[\$i]/rmsd-hist_Run\$ARGV[0].dat norm\n";
  print OUT "hist dist_O2'-O4'%\$i,0,10,0.1 out ./\$temp[\$i]/distO2O4-hist_Run\$ARGV[0].dat\n";
  print OUT "hist COM%\$i,0,10,0.1 out ./\$temp[\$i]/COM-hist_Run\$ARGV[0].dat norm\n"; 
  print OUT "hist g1%\$i,-180,180,1 out ./\$temp[\$i]/dihed_g1-hist_Run\$ARGV[0].dat norm\n";
  print OUT "hist d1%\$i,-180,180,1 out ./\$temp[\$i]/dihed_d1-hist_Run\$ARGV[0].dat norm\n";
  print OUT "hist e1%\$i,-180,180,1 out ./\$temp[\$i]/dihed_e1-hist_Run\$ARGV[0].dat norm\n";
  print OUT "hist z1%\$i,-180,180,1 out ./\$temp[\$i]/dihed_z1-hist_Run\$ARGV[0].dat norm\n";
  print OUT "hist c1%\$i,-180,180,1 out ./\$temp[\$i]/dihed_c1-hist_Run\$ARGV[0].dat norm\n";
  print OUT "hist a2%\$i,-180,180,1 out ./\$temp[\$i]/dihed_a2-hist_Run\$ARGV[0].dat norm\n";
  print OUT "hist b2%\$i,-180,180,1 out ./\$temp[\$i]/dihed_b2-hist_Run\$ARGV[0].dat norm\n";
  print OUT "hist g2%\$i,-180,180,1 out ./\$temp[\$i]/dihed_g2-hist_Run\$ARGV[0].dat norm\n";
  print OUT "hist d2%\$i,-180,180,1 out ./\$temp[\$i]/dihed_d2-hist_Run\$ARGV[0].dat norm\n";
  print OUT "hist c2%\$i,-180,180,1 out ./\$temp[\$i]/dihed_c2-hist_Run\$ARGV[0].dat norm\n";
  print OUT "hist d1%\$i,-180,180,1 COM%\$i,0,10,0.1 free \$temp[\$i] out ./\$temp[\$i]/2dhist_d1-COMene.gnu\n";
  print OUT "hist d2%\$i,-180,180,1 COM%\$i,0,10,0.1 free \$temp[\$i] out ./\$temp[\$i]/2dhist_d2-COMene.gnu\n";
  print OUT "hist c1%\$i,-180,180,1 COM%\$i,0,10,0.1 free \$temp[\$i] out ./\$temp[\$i]/2dhist_c1-COMene.gnu\n";
  print OUT "hist c2%\$i,-180,180,1 COM%\$i,0,10,0.1 free \$temp[\$i] out ./\$temp[\$i]/2dhist_c2-COMene.gnu\n";
  print OUT "hist dist_O2'-O4'%\$i,0,10,0.1 COM%\$i,0,10,0.1 free \$temp[\$i] out ./\$temp[\$i]/2dhist_d1-COMene.gnu\n";
}

system "cpptraj -i pt-ensemble.in > pt-ensemble.log";
EOF

#perl ensemble.perl $1 
cd ../ 

# stacking analysis in different temperatures
cd stack-temp
cat << EOF > analyze-stacking.perl
#!/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,C2,N3,C4,C5,C6,N6,N7,C8,N9 mindist\n";
  print OUT "distance COM :1\@N1,C6,H6,C5,H5,C4,N4,H41,H42,N3,C2,O2 :2\@N1,C2,H2,N3,C4,C5,C6,N6,H61,H62,N7,C8,H8,N9\n\n";

  print OUT "vector v_base2 :2\@N1,C2,N3,C4,C5,C6,N6,N7,C8,N9 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";

EOF

#perl analyze-stacking.perl $1 &
cd ../

# Analyzing the Suites at 298K
cd ./suites-298
cat << EOF > suites-analyzer.perl
#!/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);

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

print OUT "parm ../../build/full.topo.hmr [traj]\n\nstrip :WAT,Na+,K+,Cl-\n\nautoimage origin\n\n";

foreach \$p (1..\$ARGV[0]) {
  print OUT "trajin ../../traj.\$p.01 remdtraj remdtrajtemp 298.41 trajnames ";
  foreach \$k (1..\$#temp) {
    \$j=sprintf "%02.0f", \$k+1;
    print OUT "../../traj.\$p.\$j,";
  }
  print OUT "\n\n";
}

print OUT "trajout 298.nc netcdf\n";
open OUT2, "> pt-parmstrip.in";
print OUT2 "parm ../../build/full.topo.hmr\nparmstrip :WAT,Cl-,Na+\nparmwrite out ../../build/noWAt.topo.hmr\n";

system "cpptraj -i pt-298.in > pt-298.log";
system "cpptraj -i pt-parmstrip.in > pt-parmstrip.log";
system "perl ~/scripts/suitename/suite-classifier-niel.pl ../../build/noWAt.topo.hmr 298.nc > summary-suites.dat";
EOF


cat << EOF > pt-parmstrip.in
parm ../../build/full.topo.hmr
parmstrip :WAT,Cl-,Na+
parmwrite out ../../build/noWAt.topo.hmr
EOF

#perl suites-analyzer.perl $1
cd ../

# H-bond analysis at 298K
cd hbond-298
cat << EOF > hbond.perl
#!/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);

open OUTwat, "> pt-hb-wat.in";

print OUTwat "parm ../../build/full.topo.hmr [traj]\n\nautoimage origin\n\n";

foreach \$p (1..\$ARGV[0]) {
  print OUTwat "trajin ../../traj.\$p.01 remdtraj remdtrajtemp 298.41 trajnames ";
  foreach \$k (1..\$#temp) {
    \$j=sprintf "%02.0f", \$k+1;
    print OUTwat "../../traj.\$p.\$j,";
  }
  print OUTwat "\n\n";
}
print OUTwat "hbond nointramol solvout avghb-wat.dat solventdonor :WAT acceptormask :1-2 printatomnum dist 3 printatomnum\n";

open OUTrna, "> pt-hb-rna.in";

print OUTrna "parm ../../build/full.topo.hmr [traj]\n\nautoimage origin\n\n";

foreach \$p (1..\$ARGV[0]) {
  print OUTrna "trajin ../../traj.\$p.01 remdtraj remdtrajtemp 298.41 trajnames ";
  foreach \$k (1..\$#temp) {
    \$j=sprintf "%02.0f", \$k+1;
    print OUTrna "../../traj.\$p.\$j,";
  }
  print OUTrna "\n\n";
}
print OUTrna "autoimage origin\n";
print OUTrna "hbond :1-2 avgout avghb-rna.dat printatomnum\n";

system "cpptraj -i pt-hb-wat.in > pt-hb-wat.log";
system "cpptraj -i pt-hb-rna.in > pt-hb-rna.log";
EOF

#perl hbond.perl $1

#Analyzing cluster populations

cd ../clusters-6
cat << EOF > cluster.perl
#!/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);

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

print OUT "parm ../../build/full.topo.hmr [traj]\n\n";

foreach \$p (1..\$ARGV[0]) {
  print OUT "trajin ../../traj.\$p.01 remdtraj remdtrajtemp 298.41 trajnames ";
  foreach \$k (1..\$#temp) {
    \$j=sprintf "%02.0f", \$k+1;
    print OUT "../../traj.\$p.\$j,";
  }
  print OUT "\n\n";
}

print OUT "strip :WAT,Na+,Cl-\nrmsd :1,2@N1,C2,H2,N3,C4,C5,C6,N6,H61,H62,N7,C8,H8,N9 first mass\n\n";

print OUT "cluster summary summary.dat repout rep repfmt pdb sieve 100 clusterfmt netcdf averagelinkage rms mass :1-2 clusters 6 cpopvtime cluster_pop.agr normframe\n";

system "cpptraj.OMP -i pt-cluster.in >  pt-cluster.log";
EOF


cd ../../

cat << EOF > cpptraj_analysis-1.runit
#!/bin/tcsh
#SBATCH -p compute
#SBATCH -J $2
#SBATCH -o logtest.o%j
#SBATCH -N 1
#SBATCH --ntasks-per-node=24
#SBATCH -t 03:00:00


# Add location of Amber Python modules to default Python search path
setenv AMBERHOME "/oasis/projects/nsf/slc216/droe/GIT/amber"
set path = (\$AMBERHOME/bin \$path)

cd analysis/rep-rmsd
perl rep-rmsd.perl $1 &

cd ../remlog
perl remlog-analysis.perl $1 &

cd ../suites-298 
perl suites-analyzer.perl $1 &

cd ../hbond-298
perl hbond.perl $1 

cd ../conf-temp
perl ensemble.perl $1

cd ../stack-temp
perl analyze-stacking.perl $1

cd ../clusters-6
perl cluster.perl $1
EOF

qsub cpptraj_analysis-1.runit