#!/bin/tcsh

# first argument is the number of runs

foreach clust (0 1 2 3 4) 
#rm -r c$clust/
mkdir c$clust/
cd c$clust

cat << EOF > pt-analysis1.in
parm ../../../tip3p/ff12sb/run1/build/noWAt.topo.hmr

trajin ../ctraj.c$clust

reference ../../../A-form.pdb [a-form]

rmsd rms_2_aform :1,2&!@H= ref [a-form] mass 

distance dist_sugar_ox :1@O2' :2@O4' 
distance dist_N3-O5p :1@N3 :1@O5'
distance COM :2@N1,C2,H2,N3,C4,C5,C6,N2,H21,H22,N7,C8,H8,N9,O6 :1@N1,C6,H6,C5,H5,C4,N4,H41,H42,N3,C2,O2

dihedral g1 :1@O5' :1@C5' :1@C4' :1@C3'  type gamma
dihedral d1 :1@C5' :1@C4' :1@C3' :1@O3'  type delta
dihedral e1 :1@C4' :1@C3' :1@O3' :2@P    type epsilon
dihedral z1 :1@C3' :1@O3' :2@P   :2@O5'  type zeta
dihedral c1 :1@O4' :1@C1' :1@N1  :1@C2   type chi

dihedral a2 :1@O3' :2@P   :2@O5' :2@C5'  type alpha
dihedral b2 :2@P   :2@O5' :2@C5' :2@C4'  type beta
dihedral g2 :2@O5' :2@C5' :2@C4' :2@C3'  type gamma
dihedral d2 :2@C5' :2@C4' :2@C3' :2@O3'  type delta
dihedral c2 :2@O4' :2@C1' :2@N9  :2@C4   type chi
run

avg rms_2_aform out avg-rmsd-all.dat
hist rms_2_aform,0,5,0.1 out rmsd-hist.dat name RMSD norm
hist COM,0,10,0.1 out COM-hist.dat name CoM norm
hist dist_sugar_ox,0,10,0.1 out distO2O4-hist.dat norm
hist dist_N3-O5p,0,10,0.1 out distN3-O5p-hist.dat norm
hist rms_2_aform,0,5,0.1 dist_sugar_ox,0,10,0.1 out 2dhist_rms-dist.gnu norm
hist g1,-180,180,1 out dihed_g1-hist.dat name gamma_1 norm
hist d1,-180,180,1 out dihed_d1-hist.dat name delta_1 norm
hist e1,-180,180,1 out dihed_e1-hist.dat name epsilon_1 norm
hist z1,-180,180,1 out dihed_z1-hist.dat name zeta_1 norm
hist c1,-180,180,1 out dihed_c1-hist.dat name chi_1 norm
hist a2,-180,180,1 out dihed_a2-hist.dat name alpha_2 norm
hist b2,-180,180,1 out dihed_b2-hist.dat name beta_2 norm
hist g2,-180,180,1 out dihed_g2-hist.dat name gamma_2 norm
hist d2,-180,180,1 out dihed_d2-hist.dat name delta_2 norm
hist c2,-180,180,1 out dihed_c2-hist.dat name chi_2 norm
EOF

cat << EOF > analyze-stacking.perl
#!/usr/bin/perl

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

print OUT "parm ../../../tip3p/ff12sb/run1/build/noWAt.topo.hmr [traj]\n\nstrip :WAT,Na+,Cl-\n\nautoimage origin\n\n";
print OUT "trajin ../ctraj.c$clust\n";
print OUT "rmsd :1,2&!\@H= first\n\n";

print OUT "nativecontacts :2\@N1,N2,C2,N3,C4,C5,C6,O6,N7,C8,N9 :1\@N1,C6,C5,C4,N4,N3,C2,O2 mindist\n";
print OUT "distance COM :2\@N1,C2,H2,N3,C4,C5,C6,N2,H21,H22,N7,C8,H8,N9,O6 :1\@N1,C6,H6,C5,H5,C4,N4,H41,H42,N3,C2,O2\n\n";

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

print OUT "writedata criteria-raw.dat Contacts_00002[mindist] COM normalangle\n";

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

open INP, "< criteria-raw.dat";
open OUT2, "> stacking.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";
close INP;
close OUT2;
system "rm criteria*.dat";

EOF

# H-bond analysis at 298K
cat << EOF > hbond.perl
#!/usr/bin/perl

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

print OUTwat "parm ../../../tip3p/ff12sb/run1/build/noWAt.topo.hmr [traj]\n\nautoimage origin\n\n";
print OUTwat "trajin ../ctraj.c$clust\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 ../../../tip3p/ff12sb/run1/build/noWAt.topo.hmr [traj]\n\nautoimage origin\n\n";

print OUTrna "trajin ../ctraj.c$clust\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

perl ~/scripts/suitename/suite-classifier-niel.pl ../../../tip3p/ff12sb/run1/build/noWAt.topo.hmr ../ctraj.c$clust > summary-suites.dat

perl analyze-stacking.perl

cpptraj -i pt-analysis1.in > pt-analysis1.log
cd ../

end