#! /bin/perl
# phe_cent_3.pl Calculates the coordinates of the ring centroid and 
# 	        the distances between all PHE residues in a pdb file.
#               3/31/00 by CJ McKnight
#
# Inovke with: phe_cent_3.pl < YOUR_PDB_FILENAME
#
# As written, this script only works for pdb files with the chain specifier
#
while (<>) {
   split;

   # find the  CA atom of the phes to keep track of how many there are
   if (($_[3]=~/PHE/) && ($_[2]=~/CA/)) {
         @residue[$i] = $_[5]; 
       #if there is no chain specifer in the pdb file replace the previous line with the following line
         #@residue[$i] = $_[4];
         $i++;
   }
 
   # find the heavy atoms of the phe rings and sum the x, y and z coordinates
   if (($_[3]=~/PHE/)&&(($_[2]=~/CG/) || ($_[2]=~/CD/) || ($_[2]=~/CE/) || ($_[2]=~/CZ/))) {
         @x[$i-1] = @x[$i-1] + $_[6];
         @y[$i-1] = @y[$i-1] + $_[7];
         @z[$i-1] = @z[$i-1] + $_[8];
      #if there is no chain specifier comment out the previous 3 lines & uncomment the following 3 lines
         #@x[$i-1] = @x[$i-1] + $_[5];
         #@y[$i-1] = @y[$i-1] + $_[6];
         #@z[$i-1] = @z[$i-1] + $_[7];
  }
}

#print out the coordinates of the PHE centers
printf "\n";
for ($i = 0; $i <=  $#residue; $i++) {
   printf "coord of F";
   printf @residue[$i] ;
   printf " ring centroid = ";
   printf "%6.3f ", @x[$i]/6;
   printf "%6.3f ", @y[$i]/6;
   printf "%6.3f ", @z[$i]/6;
   printf "\n";
}

#print out the distance between the PHE residues
printf "\n";
for ($i = 0; $i <= ($#residue-1); $i++) {
   for ($j = ($i+1); $j <= ($#residue); $j++) {
      printf "The distance from the center of F";
      printf  @residue[$i] ; 
      printf " to F"; 
      printf @residue[$j] ; 
      printf " is ";
      printf "%6.3f ",(((@x[$i]/6 -@x[$j]/6)**2 + (@y[$i]/6-@y[$j]/6)**2 + (@z[$i]/6-@z[$j]/6)**2)**0.5);
      printf "Ang \n";
   }
}