#!/usr/local/bin/perl -w

# this is a filter script to move and rotate a selected set of residues.
# The residues affected are specified by the two command-line arguments.
# 	move-residues.pl start stop < decoys/T0130.try33-opt.pdb > decoys/T0130-moved.pdb
# where start is the index of the first residue to move, and stop is the
# index of the last residue to move.


# shift_vector is the translation amount by which the selected atoms are moved
# after being rotated.
@shift_vector = (-2, 0, 0);

# rotate_vector is the vector (from the center of mass) by which the selected
# points are rotated.  Note that this vector is specified in x,y,z, but
# the base of this vector is the center of mass of the selected residues
@rotate_vector = (0, 0, 1);

# theta is the angle (in radians) by which the points are rotated.
# Note: the first number is the number of degrees to rotate.  The rest
# of this line is just a conversion to radians (atan2(1,1)*4 gives pi).
$theta = 180 * atan2(1,1)*4/180;






if (@ARGV == 0)
{
  print "\nMove-residues\n\n";
  print "Usage: move-residues.pl start stop < decoys/T0130.try33-opt.pdb\n";
  print "	> decoys/T0130-moved.pdb\n\n";
  print "This is a program to move residues within a pdb file, both by\n";
  print "translation and rotation.  The first two arguments are the start\n";
  print "and stop indicies of the residues to rotate.  STDIN is assumed to\n";
  print "be the input PDB file, and STDOUT is assumed to be the output PDB\n";
  print "file.  The shift and rotation amounts are specified in variables\n";
  print "at the beginning of this program.  The axis of rotation is a\n";
  print "vector with its base at the center of mass of the selected\n";
  print "residues.\n";
  print "For the angle of rotation, positive is counter-clockwise when\n";
  print "looking down the axis towards the origin.\n";
  print "Tip: It's often easiest to rotate along one of the three\n";
  print "cardinal axes at a time, making small changes along the way,\n";
  print "rather than trying to change everything at once.\n";
  exit(0);
}


$start = shift;
$stop = shift;

while($line=<STDIN>)
{
    if (($line =~ /^HETATM/) || ($line =~ /^CONECT/))
    {
      next;
    }


    if ($line!~/ATOM/)
    {   print $line;
    	next;
    }
    $resnum=substr($line,22,4);

    if ($resnum >= $start)
    {

      #check if we've hit the last residue to rotate.  WARNING: Don't overwrite
      # this line before we print it out
      last if ($resnum > $stop);
      push @rot_list, $line;
      next;
    }
    print $line;
}


# HANDLE THE ROTATION

# get the matrix

  #unitize the rotation vector
  $sum = $rotate_vector[0] + $rotate_vector[1] + $rotate_vector[2];

&get_matrix ($rotate_vector[0]/$sum, $rotate_vector[1]/$sum, $rotate_vector[2]/$sum);

# FIND CENTER OF MASS

%atommass = GetMassHash();

for ($i=0; $i<=$#rot_list; $i++)
{

  $x= substr($rot_list[$i],30,8);    
  $y= substr($rot_list[$i],38,8);    
  $z= substr($rot_list[$i],46,8);
  $atom= substr($rot_list[$i],13,1);

  $mass = $atommass{$atom};

  $centerx += $x * $mass;
  $centery += $y * $mass;
  $centerz += $z * $mass;

  $total_mass += $mass;
}
$centerx /= $total_mass;
$centery /= $total_mass;
$centerz /= $total_mass;

print STDERR "Found Center of Mass at ($centerx, $centery, $centerz)\n";


# ROTATE AND PRINT EACH POINT

for ($i=0; $i<=$#rot_list; $i++)
{
  # get the point
  $x= substr($rot_list[$i],30,8);    
  $y= substr($rot_list[$i],38,8);    
  $z= substr($rot_list[$i],46,8);

  # subtract the center of mass from the point
  $x -= $centerx;
  $y -= $centery;
  $z -= $centerz;

  # rotate the point
  $newx = $x*$R[0][0] + $y*$R[0][1] + $z*$R[0][2];
  $newy = $x*$R[1][0] + $y*$R[1][1] + $z*$R[1][2];
  $newz = $x*$R[2][0] + $y*$R[2][1] + $z*$R[2][2];

  # add the center of mass to the point
  $newx += $centerx;
  $newy += $centery;
  $newz += $centerz;

  # print the point
  printf "%s%8.3f%8.3f%8.3f%s", 
	substr($rot_list[$i],0,30), $newx+$shift_vector[0],
	$newy+$shift_vector[1], $newz+$shift_vector[2], substr($rot_list[$i],54);
}


if ($line)
{
  print $line;

  while($line=<STDIN>)
  {
      if (($line =~ /^HETATM/) || ($line =~ /^CONECT/))
      {
        next;
      }
      print $line;
  }
}



sub get_matrix ($$$)
{
  $x = shift;
  $y = shift;
  $z = shift;

  $R[0][0] = $x*$x + cos($theta)*(1-$x*$x);

  $R[0][1] = $x*$y*(1-cos($theta)) - $z*sin($theta);

  $R[0][2] = $x*$z*(1-cos($theta)) + $y*sin($theta);

  $R[1][0] = $y*$x*(1-cos($theta)) + $z*sin($theta);

  $R[1][1] = $y*$y + cos($theta)*(1-$y*$y);

  $R[1][2] = $y*$z*(1-cos($theta)) - $x*sin($theta);

  $R[2][0] = $z*$x*(1-cos($theta)) - $y*sin($theta);

  $R[2][1] = $z*$y*(1-cos($theta)) + $x*sin($theta);

  $R[2][2] = $z*$z + cos($theta)*(1-$z*$z);

}


sub GetMassHash
{   my %mass_weights =
	(	C => '12.01',
		O => '16.00',
		N => '14.01',
		H => '1.008',
		S => '32.07',
		D => '2.0', # Deuterium (hydrogen isotope)
		A => '15.005',  #.5N + .5O -- means uncertain O/N positioning
	);
    return( %mass_weights );
}