#!/usr/bin/perl
# Copyright (c) 2011 by Global Phasing Limited
use strict;
use Math::Trig;
my $version="1.1";
print "# output from c3endo_or_c2endo_utiltor_restraints_all.pl jiffy (v$version)\n";
print "# idea is to set up BUSTER utility torsions for RNA ribose rings to\n";
print "# be pushed to canonical C3'-endo or C2'-endo conformations.\n";
print "# Ideal torsion values from http://ndbserver.rutgers.edu/archives/proj/valence/table6.html\n";
print "#\n";
print "# Help Page: https://www.globalphasing.com/buster/wiki/index.cgi?RNApuckerRestraints#pucker\n";
print "# Author: Oliver Smart\n";
print "# E-mail: buster-develop\@globalphasing.com\n";
print "# Jiffy Date: 21 October 2011\n";
print "NOTE BUSTER_REMARK3_NOTE c3endo_or_c2endo_utiltor_restraints_all.pl (v$version) used to for utility torsions to favor RNA ribose canonical C3'-endo or C2'-endo conformations.\n";

my $sigma="5.0"; # 5 degrees is loose, use 2 degrees for tight
my @utiltors_c3endo=(
"restraints to C3'-endo",
"NOTE BUSTER_UTILTOR 1   2.8 sigma  resid:C4* resid:O4* resid:C1* resid:C2*",
"NOTE BUSTER_UTILTOR 1 335.4 sigma  resid:O4* resid:C1* resid:C2* resid:C3*",
"NOTE BUSTER_UTILTOR 1  35.9 sigma  resid:C1* resid:C2* resid:C3* resid:C4*",
"NOTE BUSTER_UTILTOR 1 324.7 sigma  resid:C2* resid:C3* resid:C4* resid:O4*",
"NOTE BUSTER_UTILTOR 1  20.4 sigma  resid:C3* resid:C4* resid:O4* resid:C1*",
"NOTE BUSTER_UTILTOR 1 204.0 sigma  resid:C5* resid:C4* resid:C3* resid:C2*",
"NOTE BUSTER_UTILTOR 1  81.0 sigma  resid:C5* resid:C4* resid:C3* resid:O3*",
"NOTE BUSTER_UTILTOR 1 201.8 sigma  resid:O4* resid:C4* resid:C3* resid:O3*",
"NOTE BUSTER_UTILTOR 1 241.4 sigma  resid:C4* resid:O4* resid:C1* resid:N1", #N1 must be changed to N9 for G
"NOTE BUSTER_UTILTOR 1  44.3 sigma  resid:O3* resid:C3* resid:C2* resid:O2*");
my @utiltors_c2endo=(
"restraints to C2'-endo",
"NOTE BUSTER_UTILTOR 1 339.2 sigma  resid:C4* resid:O4* resid:C1* resid:C2*",
"NOTE BUSTER_UTILTOR 1  35.2 sigma  resid:O4* resid:C1* resid:C2* resid:C3*",
"NOTE BUSTER_UTILTOR 1 324.6 sigma  resid:C1* resid:C2* resid:C3* resid:C4*",
"NOTE BUSTER_UTILTOR 1  24.2 sigma  resid:C2* resid:C3* resid:C4* resid:O4*",
"NOTE BUSTER_UTILTOR 1 357.7 sigma  resid:C3* resid:C4* resid:O4* resid:C1*",
"NOTE BUSTER_UTILTOR 1 263.4 sigma  resid:C5* resid:C4* resid:C3* resid:C2*",
"NOTE BUSTER_UTILTOR 1 147.3 sigma  resid:C5* resid:C4* resid:C3* resid:O3*",
"NOTE BUSTER_UTILTOR 1 268.1 sigma  resid:O4* resid:C4* resid:C3* resid:O3*",
"NOTE BUSTER_UTILTOR 1 216.6 sigma  resid:C4* resid:O4* resid:C1* resid:N1", #N1 must be changed to N9 for G
"NOTE BUSTER_UTILTOR 1 319.7 sigma  resid:O3* resid:C3* resid:C2* resid:O2*");
my $c1_match=0;   my $c2_match=0;   my $c3_match=0;   my $c4_match=0;
my $c1_resnum=""; my $c2_resnum=""; my $c3_resnum=""; my $c4_resnum="";
my $c1_chsnam=""; my $c2_chsnam=""; my $c3_chsnam=""; my $c4_chsnam="";
my @c1_coords;    my @c2_coords;    my @c3_coords;    my @c4_coords;   

while (<>) { 
  my $line = $_;
  if ($line =~/^ATOM |^HETATM/) { # have an atom line
    my $atnam =  substr($line,13,3);
    my $resnam =  substr($line,17,3);
    my $chsnam =  substr($line,21,1);
    my $resnum =  substr($line,22,4);
    my $xcoord = substr($line,30,8);
    my $ycoord = substr($line,38,8);
    my $zcoord = substr($line,46,8);
    $resnam =~ s/ //g;
    $resnum =~ s/ //g;
    if ($resnam=~/A|G|C|U/) { # have an RNA residue
      if ($atnam =~/C1\*/) { # have atom C1* store its info
         #print "#debug  C1* atnam=\"$atnam\" resnam=\"$resnam\" chsnam=\"$chsnam\" resnum=\"$resnum\" \n";
         #print "#debug  C1* coords \"$xcoord\" \"$ycoord\" \"$zcoord\"\n";
         $c1_match = 1;
         $c1_resnum = $resnum;
         $c1_chsnam = $chsnam;
         $c1_coords[0] = $xcoord;         
         $c1_coords[1] = $ycoord;         
         $c1_coords[2] = $zcoord;         
      }
      if ($atnam =~/C2\*/) { # have atom C2*
         #print "#debug  C2* atnam=\"$atnam\" resnam=\"$resnam\" chsnam=\"$chsnam\" resnum=\"$resnum\" \n";
         #print "#debug  C2* coords \"$xcoord\" \"$ycoord\" \"$zcoord\"\n";
         $c2_match = 1;
         $c2_resnum = $resnum;
         $c2_chsnam = $chsnam;
         $c2_coords[0] = $xcoord;         
         $c2_coords[1] = $ycoord;         
         $c2_coords[2] = $zcoord;         
      }
      if ($atnam =~/C3\*/) { # have atom C3*
         #print "#debug  C3* atnam=\"$atnam\" resnam=\"$resnam\" chsnam=\"$chsnam\" resnum=\"$resnum\" \n";
         #print "#debug  C3* coords \"$xcoord\" \"$ycoord\" \"$zcoord\"\n";
         $c3_match = 1;
         $c3_resnum = $resnum;
         $c3_chsnam = $chsnam;
         $c3_coords[0] = $xcoord;         
         $c3_coords[1] = $ycoord;         
         $c3_coords[2] = $zcoord;         
      }
      if ($atnam =~/C4\*/) { # have atom C4*
         #print "#debug  C4* atnam=\"$atnam\" resnam=\"$resnam\" chsnam=\"$chsnam\" resnum=\"$resnum\" \n";
         #print "#debug  C4* coords \"$xcoord\" \"$ycoord\" \"$zcoord\"\n";
         $c4_match = 1;
         $c4_resnum = $resnum;
         $c4_chsnam = $chsnam;
         $c4_coords[0] = $xcoord;         
         $c4_coords[1] = $ycoord;         
         $c4_coords[2] = $zcoord;         
      }
      if (($c1_match) && ($c2_match) && ($c3_match)  && ($c4_match)) { # have we match C1*, C2*, C3* and C4*?
          #print "# debug all 4 match\n";
          $c1_match = 0; 
          $c2_match = 0; 
          $c3_match = 0; 
          $c4_match = 0;  # clear matches
          # check that the resnum chsnam are the same
          die "ERROR have an incomplete residue at or before $chsnam $resnum" 
          if ( ($c1_resnum != $c2_resnum) ||  ($c1_resnum != $c3_resnum) || ($c1_resnum != $c4_resnum));
          die "ERROR have an incomplete residue at or before $chsnam $resnum" 
          if ( ($c1_chsnam != $c2_chsnam) ||  ($c1_chsnam != $c3_chsnam) || ($c1_chsnam != $c4_chsnam));
          # so now should have the coordinates of all C1*, C2*, C3*, C4*
          #print "# debug resnam=\"$resnam\" chsnam=\"$chsnam\" resnum=\"$resnum\"\n";
	  #print "# debug  C1* coords @c1_coords \n";
          #print "# debug  C2* coords @c2_coords \n";
          #print "# debug  C3* coords @c3_coords \n";
          #print "# debug  C4* coords @c4_coords \n";
          # find the torsion angle C1*-C2*-C3*-C4*
          my $thistord = torsion_angle_in_degs(  \@c1_coords, \@c2_coords,  \@c3_coords, \@c4_coords);
          print "\n";
          print "# resnam=\"$resnam\" chsnam=\"$chsnam\" resnum=\"$resnum\" C1*-C2*-C3*-C4* tors= ";
          printf ("%7.2f degs\n", $thistord);

          # if the torsion is positive write out C3'-endo restraints for the residue
          if ($thistord >= 0.) {
             print "# resnam=\"$resnam\" chsnam=\"$chsnam\" resnum=\"$resnum\" +ve torsion, ";
             foreach (@utiltors_c3endo) {
                my $thistors= $_;
                $thistors =~ s/sigma/$sigma/g;
                $thistors =~ s/resid/$chsnam\|$resnum/g;
                $thistors =~ s/N1/N9/ if ($resnam=~/G|A/);
                print "$thistors\n";
             }
           }
           else {
             print "# resnam=\"$resnam\" chsnam=\"$chsnam\" resnum=\"$resnum\" -ve torsion, ";
             foreach (@utiltors_c2endo) {
                my $thistors= $_;
                $thistors =~ s/sigma/$sigma/g;
                $thistors =~ s/resid/$chsnam\|$resnum/g;
                $thistors =~ s/N1/N9/ if ($resnam=~/G|A/);
                print "$thistors\n";
             }
           }
      } # have we match C1*, C2*, C3* and C4*?

    } # have an RNA residue

    #if (($atnam =~/C1\*/) && ($resnam=~/A|G|C|U/)) { # use atom C1*
    #  print "\n#  resnam=\"$resnam\" chsnam=\"$chsnam\" resnum=\"$resnum\"\n";
    #  foreach (@utiltors_c3endo) {
    #     my $thistors= $_;
    #     $thistors =~ s/sigma/$sigma/g;
    #     $thistors =~ s/resid/$chsnam\|$resnum/g;
    #     $thistors =~ s/N1/N9/ if ($resnam=~/G|A/);
    #     print "$thistors\n";
    #  }
    #  foreach (@utiltors_c2endo) {
    #     my $thistors= $_;
    #     $thistors =~ s/sigma/$sigma/g;
    #     $thistors =~ s/resid/$chsnam\|$resnum/g;
    #     $thistors =~ s/N1/N9/ if ($resnam=~/G|A/);
    #     print "$thistors\n";
    #  }
    #}

  }
  else {
    #print $line;
  }

}
sub torsion_angle_in_degs{  
   # torsion angle routine based on very old fortran tic code oss phd work
   my ($ref1,$ref2,$ref3,$ref4) = @_;
   my @xyz1 = @{$ref1};
   my @xyz2 = @{$ref2};
   my @xyz3 = @{$ref3};
   my @xyz4 = @{$ref4};
   #print "#debug call to torsion_angle_in_degs\n";
   #print "#debug xyz1=@xyz1\n";
   #print "#debug xyz2=@xyz2\n";
   #print "#debug xyz3=@xyz3\n";
   #print "#debug xyz4=@xyz4\n";

   my @veca; # vector from atom two to atom one
   $veca[0] = $xyz1[0] - $xyz2[0];
   $veca[1] = $xyz1[1] - $xyz2[1];
   $veca[2] = $xyz1[2] - $xyz2[2];
  
   my @vecb; # vector from atom two to atom three
   $vecb[0] = $xyz3[0] - $xyz2[0];
   $vecb[1] = $xyz3[1] - $xyz2[1];
   $vecb[2] = $xyz3[2] - $xyz2[2];

   my @vecc; # vector from atom two to atom four  
   $vecc[0] = $xyz4[0] - $xyz2[0];
   $vecc[1] = $xyz4[1] - $xyz2[1];
   $vecc[2] = $xyz4[2] - $xyz2[2];
 
   # u = a cross b
   my @vecu;
   $vecu[0] = $veca[1]*$vecb[2] - $veca[2]*$vecb[1];
   $vecu[1] = $veca[2]*$vecb[0] - $veca[0]*$vecb[2];
   $vecu[2] = $veca[0]*$vecb[1] - $veca[1]*$vecb[0];
   # v = c cross b
   my @vecv;
   $vecv[0] = $vecc[1]*$vecb[2] - $vecc[2]*$vecb[1];
   $vecv[1] = $vecc[2]*$vecb[0] - $vecc[0]*$vecb[2];
   $vecv[2] = $vecc[0]*$vecb[1] - $vecc[1]*$vecb[0];
   
   # unit vector u
   my $lenu = $vecu[0]*$vecu[0] + $vecu[1]*$vecu[1] + $vecu[2]*$vecu[2];
   $lenu = sqrt($lenu);
   $vecu[0] =  $vecu[0]/$lenu;
   $vecu[1] =  $vecu[1]/$lenu;
   $vecu[2] =  $vecu[2]/$lenu;

   # unit vector v
   my $lenv = $vecv[0]*$vecv[0] + $vecv[1]*$vecv[1] + $vecv[2]*$vecv[2];
   $lenv = sqrt($lenv);
   $vecv[0] =  $vecv[0]/$lenv;
   $vecv[1] =  $vecv[1]/$lenv;
   $vecv[2] =  $vecv[2]/$lenv;

   # cosphi = u dot v
   my $cosphi = $vecu[0]*$vecv[0] + $vecu[1]*$vecv[1]+ $vecu[2]*$vecv[2];

   # find cross product ucrossv = vecu cross vecv
   my @ucrossv;
   $ucrossv[0] = $vecu[1]*$vecv[2] - $vecu[2]*$vecv[1];
   $ucrossv[1] = $vecu[2]*$vecv[0] - $vecu[0]*$vecv[2];
   $ucrossv[2] = $vecu[0]*$vecv[1] - $vecu[1]*$vecv[0];
 
   # now unit vector vecb
   my $lenb = $vecb[0]*$vecb[0] + $vecb[1]*$vecb[1] + $vecb[2]*$vecb[2];
   $lenb = sqrt($lenb);
   $vecb[0] =  $vecb[0]/$lenb;
   $vecb[1] =  $vecb[1]/$lenb;
   $vecb[2] =  $vecb[2]/$lenb;

   # sinph = vecb dot ucrossv
   my $sinphi = $vecb[0]*$ucrossv[0] + $vecb[1]*$ucrossv[1] + $vecb[2]*$ucrossv[2];
   

   my $phi= atan2($sinphi,$cosphi); # will be returned in radians
   $phi = (180./pi)*$phi; # convert to degs

   #print "#debug phi=$phi     cosphi=$cosphi sinphi=$sinphi\n";
   return $phi;
}
