#!/usr/sbin/perl
#
# WARNING: THIS PROGRAM IS IN THE TESTING STAGE
# 
# program to create AMBER disang restraints for DNA 2'-ENDO ribose rings
#
# usage   make_ribose_dihedrals.pl res1 res2 res3 res4 outfile
#
# EXAMPLE for a 12-mer duplex: 
#
#	./make_ribose_dihedrals.pl 1 12 13 24 myBdna.disang	
#
#
$start1 = shift(@ARGV);
$stop1 = shift(@ARGV);
$start2 = shift(@ARGV);
$stop2 = shift(@ARGV);
$outfile = shift(@ARGV);
#
#
#	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";
#
#
for $count ($start1..$stop1)
  {
  $previous = $count - 1;
  $next = $count + 1;
  print OUT "\n";
  print OUT "####################################################### \n";
  print OUT "#   B-FORM DEOXY-RIBOSE RESTRAINTS FOR RESIDUE $count \n";
  print OUT "####################################################### \n";
  print OUT "\n";
  if ($count - $start1 != 0)    # DONT DO ALPHA,BETA FOR FIRST RESIDUE
    {
    print OUT "# Torsional Angle: alpha = -46 \n";
    print OUT " &rst \n";
    print OUT " iat(1) =   $previous , iat(2) =   $count , iat(3) =   $count , iat(4) =   $count , \n";
    print OUT " atnam(1) = 'O3''' , atnam(2) = 'P   ', atnam(3) = 'O5'''  , \n";
    print OUT " atnam(4) = 'C5''' , iresid = 1, \n";
    print OUT " r1  =  -56.0 , r2  =  -51.0 , r3  =  -41.0 , r4  =  -31.0 , \n";
    print OUT " rk2 =  500.0 , rk3 =  500.0 \n";
    print OUT " &end \n";
    print OUT "# Torsional Angle: beta = -147 \n";
    print OUT " &rst \n";
    print OUT " iat(1) =   $count , iat(2) =   $count , iat(3) =   $count , iat(4) =   $count , \n";
    print OUT " atnam(1) = 'P   ' , atnam(2) = 'O5''', atnam(3) = 'C5'''  , \n";
    print OUT " atnam(4) = 'C4''' , iresid = 1, \n";
    print OUT " r1  = -157.0 , r2  = -152.0 , r3  = -142.0 , r4  = -137.0 , \n";
    print OUT " rk2 =  500.0 , rk3 =  500.0 \n";
    print OUT " &end \n";
    }
#
#    DO GAMMA,DELTA FOR ALL BASES OF STRAND 1
#
  print OUT "\n";
  print OUT "# Torsional Angle: gamma = 36 \n";
  print OUT " &rst \n";
  print OUT " iat(1) =   $count , iat(2) =   $count , iat(3) =   $count , iat(4) =   $count , \n";
  print OUT " atnam(1) = 'O5''' , atnam(2) = 'C5''', atnam(3) = 'C4'''  , \n";
  print OUT " atnam(4) = 'C3''' , iresid = 1, \n";
  print OUT " r1  =   26.0 , r2  =   31.0 , r3  =   41.0 , r4  =   46.0 , \n";
  print OUT " rk2 =  500.0 , rk3 =  500.0 \n";
  print OUT " &end \n";
  print OUT "# Torsional Angle: delta = 157 \n";
  print OUT " &rst \n";
  print OUT " iat(1) =   $count , iat(2) =   $count , iat(3) =   $count , iat(4) =   $count , \n";
  print OUT " atnam(1) = 'C5''' , atnam(2) = 'C4''', atnam(3) = 'C3'''  , \n";
  print OUT " atnam(4) = 'O3''' , iresid = 1, \n";
  print OUT " r1  =  147.0 , r2  =  152.0 , r3  =  162.0 , r4  =  167.0 , \n";
  print OUT " rk2 =  500.0 , rk3 =  500.0 \n";
  print OUT " &end \n";
#
#   DO EPSILON AND ZETA EXCEPT FOR THE LAST RIBOSE OF FIRST STRAND
#
  if ($stop1 - $count != 0)
  {
  print OUT "# Torsional Angle: epsilon = 155 \n";
  print OUT " &rst \n";
  print OUT " iat(1) =   $count , iat(2) =   $count , iat(3) =   $count , iat(4) =   $next , \n";
  print OUT " atnam(1) = 'C4''' , atnam(2) = 'C3''', atnam(3) = 'O3'''  , \n";
  print OUT " atnam(4) = 'P   ' , iresid = 1, \n";
  print OUT " r1  =  145.0 , r2  =  150.0 , r3  =  160.0 , r4  =  165.0 , \n";
  print OUT " rk2 =  500.0 , rk3 =  500.0 \n";
  print OUT " &end \n";
  print OUT "# Torsional Angle: zeta = -96 \n";
  print OUT " &rst \n";
  print OUT " iat(1) =   $count , iat(2) =   $count , iat(3) =   $next , iat(4) =   $next , \n";
  print OUT " atnam(1) = 'C3''' , atnam(2) = 'O3''', atnam(3) = 'P   '  , \n";
  print OUT " atnam(4) = 'O5''' , iresid = 1, \n";
  print OUT " r1  = -106.0 , r2  = -101.0 , r3  =  -91.0 , r4  =  -86.0 , \n";
  print OUT " rk2 =  500.0 , rk3 =  500.0 \n";
  print OUT " &end \n";
  }
}
# END OF FIRST STRAND
###############################
# REPEAT FOR SECOND STRAND
###############################
for $count ($start2..$stop2)
  {
  $previous = $count - 1;
  $next = $count + 1;
  print OUT "\n";
  print OUT "#   B-FORM DEOXY-RIBOSE RESTRAINTS FOR RESIDUE $count \n";
  print OUT "\n";
  if ($count - $start2 != 0)    # DONT DO ALPHA,BETA FOR FIRST RESIDUE OF A STRAND
    {
    print OUT "# Torsional Angle: alpha = -46 \n";
    print OUT " &rst \n";
    print OUT " iat(1) =   $previous , iat(2) =   $count , iat(3) =   $count , iat(4) =   $count , \n";
    print OUT " atnam(1) = 'O3''' , atnam(2) = 'P   ', atnam(3) = 'O5'''  , \n";
    print OUT " atnam(4) = 'C5''' , iresid = 1, \n";
    print OUT " r1  =  -56.0 , r2  =  -51.0 , r3  =  -41.0 , r4  =  -36.0 , \n";
    print OUT " rk2 =  500.0 , rk3 =  500.0 \n";
    print OUT " &end \n";
    print OUT "# Torsional Angle: beta = -147 \n";
    print OUT " &rst \n";
    print OUT " iat(1) =   $count , iat(2) =   $count , iat(3) =   $count , iat(4) =   $count , \n";
    print OUT " atnam(1) = 'P   ' , atnam(2) = 'O5''', atnam(3) = 'C5'''  , \n";
    print OUT " atnam(4) = 'C4''' , iresid = 1, \n";
    print OUT " r1  = -157.0 , r2  = -152.0 , r3  = -142.0 , r4  = -137.0 , \n";
    print OUT " rk2 =  500.0 , rk3 =  500.0 \n";
    print OUT " &end \n";
    }
#
#    DO GAMMA,DELTA FOR ALL BASES OF STRAND 2
#
  print OUT "\n";
  print OUT "# Torsional Angle: gamma = 36 \n";
  print OUT " &rst \n";
  print OUT " iat(1) =   $count , iat(2) =   $count , iat(3) =   $count , iat(4) =   $count , \n";
  print OUT " atnam(1) = 'O5''' , atnam(2) = 'C5''', atnam(3) = 'C4'''  , \n";
  print OUT " atnam(4) = 'C3''' , iresid = 1, \n";
  print OUT " r1  =   26.0 , r2  =   31.0 , r3  =   41.0 , r4  =   46.0 , \n";
  print OUT " rk2 =  500.0 , rk3 =  500.0 \n";
  print OUT " &end \n";
  print OUT "# Torsional Angle: delta = 157 \n";
  print OUT " &rst \n";
  print OUT " iat(1) =   $count , iat(2) =   $count , iat(3) =   $count , iat(4) =   $count , \n";
  print OUT " atnam(1) = 'C5''' , atnam(2) = 'C4''', atnam(3) = 'C3'''  , \n";
  print OUT " atnam(4) = 'O3''' , iresid = 1, \n";
  print OUT " r1  =  147.0 , r2  =  152.0 , r3  =  162.0 , r4  =  167.0 , \n";
  print OUT " rk2 =  500.0 , rk3 =  500.0 \n";
  print OUT " &end \n";
#
#   DO EPSILON AND ZETA EXCEPT FOR THE LAST RIBOSE OF SECOND STRAND
#
  if ($stop2 - $count != 0)
  {
  print OUT "# Torsional Angle: epsilon = 155 \n";
  print OUT " &rst \n";
  print OUT " iat(1) =   $count , iat(2) =   $count , iat(3) =   $count , iat(4) =   $next , \n";
  print OUT " atnam(1) = 'C4''' , atnam(2) = 'C3''', atnam(3) = 'O3'''  , \n";
  print OUT " atnam(4) = 'P   ' , iresid = 1, \n";
  print OUT " r1  =  145.0 , r2  =  150.0 , r3  =  160.0 , r4  =  165.0 , \n";
  print OUT " rk2 =  500.0 , rk3 =  500.0 \n";
  print OUT " &end \n";
  print OUT "# Torsional Angle: zeta = -96 \n";
  print OUT " &rst \n";
  print OUT " iat(1) =   $count , iat(2) =   $count , iat(3) =   $next , iat(4) =   $next , \n";
  print OUT " atnam(1) = 'C3''' , atnam(2) = 'O3''', atnam(3) = 'P   '  , \n";
  print OUT " atnam(4) = 'O5''' , iresid = 1, \n";
  print OUT " r1  = -106.0 , r2  = -101.0 , r3  =  -91.0 , r4  =  -68.0 , \n";
  print OUT " rk2 =  500.0 , rk3 =  500.0 \n";
  print OUT " &end \n";
  }
}
# END OF SECOND STRAND

#  print OUT "\n &rst iat(1) = 0 &end\n";