#!/usr/sbin/perl
# 
# program to create AMBER disang restraints for alpha helices and beta sheets
#
# D. E. Volk, University of Texas Medical Branch, Galveston, TX.
#  
#
# usage   make_protein_dihedrals.pl infile outfile
#
#	EXAMPLE infile contains only these lines
#	helix 4 10
#	sheet 15   21
#	sheet 27 33
#
#
$infile = shift(@ARGV);
$outfile = shift(@ARGV);
#
#	OPEN INPUT FILE
#
open(INFILE, "$infile") || die "Cannot open file $infile";
#
#	REMOVE OLD OUTPUT FILE AND CREATE NEW ONE
#
if (-e "$outfile") {
   unlink("$outfile");
   }
open(OUT, ">$outfile") || die "Cannot open file $outfile";
#
#	PRINT TWO LINE HEADER OF DISANG FILE
#
print OUT " SANDER DIHEDRAL RESTRAINT INPUT FILE\n";
print OUT "      ******** Section Three ********\n\n";
#
#	READ infile AND WRITE AMBER RESTRAINTS TO outfile
#
while (<INFILE>)
  {
   chomp;    #get rid of newline
   s/^\s+//;
   ($shape,$first,$second)=split(/\s+/,$_);
   if ($shape eq 'helix') 
     {
        print OUT "#   ALPHA HELIX PHI RESTRAINTS FOR RESIDUES $first - $second \n";
      for $count ($first..$second)
       {
        $res1 = ($count-1);
        $res2 = $count;
        print OUT "\n";
        print OUT "# helix phi for residue $count\n";
	print OUT " &rst\n";
	print OUT " iat(1) =  $res1     , iat(2) =  $res2     ,\n";
	print OUT " iat(3) =  $res2     , iat(4) =  $res2     ,\n";
	print OUT " atnam(1) = 'C   ' , atnam(2) = 'N   ' , atnam(3) = 'CA  ' ,\n";
        print OUT " atnam(4) = 'C   ' , iresid = 1,\n";
#	print OUT " r1  = -110.0 , r2  =  -90.0 , r3  =  -30.0 , r4  =  -10.0,\n";
	print OUT " r1  =  -80.0 , r2  =  -70.0 , r3  =  -50.0 , r4  =  -40.0,\n";
	print OUT " rk2 =  500.0 , rk3 =  500.0\n"; 
        print OUT " &end\n";
       }
	print OUT "\n";
        print OUT "#    ALPHA HELIX PSI RESTRAINTS FOR RESIDUES $first - $second \n";
      for $count ($first..$second)
       {
        $res1 = $count;
        $res2 = ($count+1);
        print OUT "\n";
        print OUT "# helix psi for residue $count\n";
        print OUT " &rst\n";
        print OUT " iat(1) =  $res1     , iat(2) =  $res1     ,\n";
        print OUT " iat(3) =  $res1     , iat(4) =  $res2     ,\n";
        print OUT " atnam(1) = 'N   ' , atnam(2) = 'CA  ' , atnam(3) = 'C   ' ,\n";
        print OUT " atnam(4) = 'N   ' , iresid = 1,\n";
#       print OUT " r1  =  -90.0 , r2  =  -70.0 , r3  =  -10.0 , r4  =   10.0,\n";
        print OUT " r1  =  -60.0 , r2  =  -50.0 , r3  =  -30.0 , r4  =  -20.0,\n";
        print OUT " rk2 =  500.0 , rk3 =  500.0\n";
        print OUT " &end\n";
       }
        print OUT "\n";
        print OUT "#    ALPHA HELIX H-BOND RESTRAINTS FOR RESIDUES $first - $second \n";
      $stop = ($second-4);
      for $count ($first..$stop)
       {
        $res1 = $count;
        $res2 = ($count+4);
        print OUT "\n";
        print OUT "# helical h-bond restraint for residue $count\n";
        print OUT " &rst\n";
        print OUT " iat(1) =  $res1     , iat(2) =  $res2     ,\n";
        print OUT " atnam(1) = 'O   ' , atnam(2) = 'H   ' ,\n";
        print OUT " iresid = 1,\n";
        print OUT " r1  =  1.75 , r2  =   1.80 , r3  =  2.20  , r4  =   2.40,\n";
        print OUT " rk2 =  500.0 , rk3 =  500.0\n";
        print OUT " &end\n";
       }


     }
     elsif ($shape eq 'sheet')
     {
        print OUT "#   BETA SHEET PHI RESTRAINTS FOR RESIDUES $first - $second \n";
      for $count ($first..$second)
       {
        $res1 = ($count-1);
        $res2 = $count;
        print OUT "\n";
        print OUT "# beta-sheet phi restraint for residue $count\n";
        print OUT " &rst\n";
        print OUT " iat(1) =  $res1     , iat(2) =  $res2     ,\n";
        print OUT " iat(3) =  $res2     , iat(4) =  $res2     ,\n";
        print OUT " atnam(1) = 'C   ' , atnam(2) = 'N   ' , atnam(3) = 'CA  ' ,\n";
        print OUT " atnam(4) = 'C   ' , iresid = 1,\n";
#       print OUT " r1  = -160.0 , r2  = -140.0 , r3  =  -80.0 , r4  =  -60.0,\n";
        print OUT " r1  = -130.0 , r2  = -120.0 , r3  = -100.0 , r4  =  -90.0,\n";
        print OUT " rk2 =  500.0 , rk3 =  500.0\n";
        print OUT " &end\n";
       }
        print OUT "\n";
        print OUT "#    BETA SHEET PSI RESTRAINTS FOR RESIDUES $first - $second \n";
      for $count ($first..$second)
       {
        $res1 = $count;
        $res2 = ($count+1);
        print OUT "\n";
        print OUT "# beta-sheet psi restraint for residue $count\n";
        print OUT " &rst\n";
        print OUT " iat(1) =  $res1     , iat(2) =  $res1     ,\n";
        print OUT " iat(3) =  $res1     , iat(4) =  $res2     ,\n";
        print OUT " atnam(1) = 'N   ' , atnam(2) = 'CA  ' , atnam(3) = 'C   ' ,\n";
        print OUT " atnam(4) = 'N   ' , iresid = 1,\n";
#       print OUT " r1  =   90.0 , r2  =  110.0 , r3  =  170.0 , r4  =  190.0,\n";
        print OUT " r1  =  120.0 , r2  =  130.0 , r3  =  150.0 , r4  =  160.0,\n";
        print OUT " rk2 =  500.0 , rk3 =  500.0\n";
        print OUT " &end\n";
       }
     }
    else {die "$shape is not helix or sheet".\n};
   }
   print OUT "\n &rst iat(1) = 0 &end\n";