#!/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\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 "strip :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,O4,N3,H3,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"; 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,O4,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,O4,N3,H3,C2,O2 :2\@N1,C2,H2,N3,C4,C5,C6,N6,H61,H62,N7,C8,H8,N9\n\n"; print OUT "vector v_base1 :1\@N1,C6,C5,C4,O4,N3,C2,O2 corrplane\n"; print OUT "vector v_base2 :2\@N1,C2,N3,C4,C5,C6,N6,N7,C8,N9 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_00002[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 () { @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 cd ../clusters-6 #Analyzing cluster populations 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 -A cheatham-lp #SBATCH -p cheatham-lp #SBATCH -J $2 #SBATCH -o logtest.o%j #SBATCH -N 1 #SBATCH -n 8 #SBATCH -t 10:00:00 setenv PATH "/uufs/chpc.utah.edu/sys/installdir/openmpi/1.8.4i/bin:\$PATH" if ( -z "\${LD_LIBRARY_PATH}" ) then setenv LD_LIBRARY_PATH "/uufs/chpc.utah.edu/sys/installdir/openmpi/1.8.4i/lib" else setenv LD_LIBRARY_PATH "/uufs/chpc.utah.edu/sys/installdir/openmpi/1.8.4i/lib:\$LD_LIBRARY_PATH" endif source /uufs/chpc.utah.edu/common/home/u0827715/Amber/GIT/amber-lonepeak/amber.csh setenv AMBERHOME /uufs/chpc.utah.edu/common/home/u0827715/Amber/GIT/amber-lonepeak 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