#!/usr/bin/perl -w

# Calculate *.sc.pdb files from *.pdb files
# type ./calSc.pl -if filename to run

if(!@ARGV){
	printf "calSc.pl: generate pdb files with side-chain centers from pdb files with only backbone atoms (including CB).\n";
	printf "example: calSc.pl -if 1a32.pdb\n";
	exit();
}
$isEst=0;   # if ==1: estimate the sc position from C_i-1, C_i, and C_i+1
$pi=3.1415926;
# get command line arguments
foreach $i (0..$#ARGV) {
    if($ARGV[$i] eq "-if") {
          $InFile=$ARGV[$i+1];
          $i++;
    }
    if($ARGV[$i] eq "-of") {
        open OUT, ">$ARGV[$i+1]";
        $i++;
    }
    if($ARGV[$i] eq "-e") {
          $isEst=1;
    }
}

@R_sc=(1.5,2.0,2.5,3.1,3.4,0.0,3.1,2.3,3.5,2.6,3.0,2.5,1.9,3.1,4.1,1.9
,1.9,2.0,3.9,3.8);

$ResLab=0;
$Res{"ALA"}=$ResLab;$Res{$ResLab++}="ALA";$Res{"CYS"}=$ResLab;$Res{$ResLab++}="CYS";
$Res{"ASP"}=$ResLab;$Res{$ResLab++}="ASP";$Res{"GLU"}=$ResLab;$Res{$ResLab++}="GLU";
$Res{"PHE"}=$ResLab;$Res{$ResLab++}="PHE";$Res{"GLY"}=$ResLab;$Res{$ResLab++}="GLY";
$Res{"HIS"}=$ResLab;$Res{$ResLab++}="HIS";$Res{"ILE"}=$ResLab;$Res{$ResLab++}="ILE";
$Res{"LYS"}=$ResLab;$Res{$ResLab++}="LYS";$Res{"LEU"}=$ResLab;$Res{$ResLab++}="LEU";
$Res{"MET"}=$ResLab;$Res{$ResLab++}="MET";$Res{"ASN"}=$ResLab;$Res{$ResLab++}="ASN";
$Res{"PRO"}=$ResLab;$Res{$ResLab++}="PRO";$Res{"GLN"}=$ResLab;$Res{$ResLab++}="GLN";
$Res{"ARG"}=$ResLab;$Res{$ResLab++}="ARG";$Res{"SER"}=$ResLab;$Res{$ResLab++}="SER";
$Res{"THR"}=$ResLab;$Res{$ResLab++}="THR";$Res{"VAL"}=$ResLab;$Res{$ResLab++}="VAL";
$Res{"TRP"}=$ResLab;$Res{$ResLab++}="TRP";$Res{"TYR"}=$ResLab;$Res{$ResLab++}="TYR";

if($InFile =~ /(.+)\.pdb$/){
$filename = $1;
}

else{$filename=$InFile;}

if($isEst==0){
    if(-e "$filename.pdb"){
	open PDBFILE, "$filename.pdb";
    }
    elsif(-e "$filename.ab.pdb"){
	open PDBFILE, "$filename.ab.pdb";
    }
    else {
	printf "Pdb file for $filename cannot not be found!\n";
	next;
    }
}
else {
    if(-e "$filename.pdb"){
	open PDBFILE, "$filename.pdb";
    }
    elsif(-e "$filename.ca.pdb"){
	open PDBFILE, "$filename.ca.pdb";
    }
    else {
	printf "Pdb file for $filename cannot not be found!\n";
	next;
    }
}
    
$lineCnt=0;
while(<PDBFILE>){
	$line[$lineCnt++]=$_;
}

close PDBFILE;

if($isEst==0){
    open OUT, ">$filename.sc.pdb";
}
else {
    open OUT, ">$filename.sce.pdb";
}
if($isEst==0){
    $ir=1;
    for($i=0;$i<$lineCnt;++$i) {
	if(substr($line[$i],16,4) =~ /[BCDE234][A-Z][A-Z][A-Z]/) {next;}
	$atomName= substr($line[$i], 12,4); $atomName=~ s/^\s+|\s+$//g; 
	$resName = substr($line[$i], 17,3); $resName =~ s/^\s+|\s+$//g;
	$resNum  = substr($line[$i], 22,4); $resNum  =~ s/^\s+|\s+$//g;
	$chainID = substr($line[$i], 21,1); $chainID = "*" if($chainID eq " ");
	$x       = substr($line[$i], 30,8); $x       =~ s/^\s+|\s+$//g;
	$y       = substr($line[$i], 38,8); $y       =~ s/^\s+|\s+$//g;
	$z       = substr($line[$i], 46,8); $z       =~ s/^\s+|\s+$//g;
	if($atomName ne "CA" && $atomName ne "CB"){next;}
	if($atomName eq "CA") {
	    printf OUT "ATOM  %5d  CA  %s %s%4d    %8.3f%8.3f%8.3f\n",$i,$resName,$chainID,$ir,$x,$y,$z;
	    if($resName eq "GLY"){$ir++;next;}
	    if($i+1!=$lineCnt) {
		$atomName2= substr($line[$i+1], 12,4); $atomName2=~ s/^\s+|\s+$//g;
		$resName2 = substr($line[$i+1], 17,3); $resName2 =~ s/^\s+|\s+$//g;
		$resNum2  = substr($line[$i+1], 22,4); $resNum2  =~ s/^\s+|\s+$//g;
		$chainID2 = substr($line[$i+1], 21,1); $chainID2 = "*" if($chainID eq " ");
		$x2  = substr($line[$i+1], 30,8); $x2  =~ s/^\s+|\s+$//g;
		$y2  = substr($line[$i+1], 38,8); $y2  =~ s/^\s+|\s+$//g;
		$z2  = substr($line[$i+1], 46,8); $z2  =~ s/^\s+|\s+$//g;
		if($atomName2 eq "CB" && $resNum2==$resNum) {
		    $dis=sqrt(($x-$x2)*($x-$x2)+($y-$y2)*($y-$y2)+($z-$z2)*($z-$z2));
		    $x3 = $R_sc[$Res{$resName}]*($x2-$x)/$dis + $x;
		    $y3 = $R_sc[$Res{$resName}]*($y2-$y)/$dis + $y;
		    $z3 = $R_sc[$Res{$resName}]*($z2-$z)/$dis + $z;
		    printf OUT "ATOM  %5d  SC  %s %s%4d    %8.3f%8.3f%8.3f\n",$i+1,$resName,$chainID,$ir,$x3,$y3,$z3;
		    ++$i;
		}
	    }
	}
	$ir++;
    }
}
else {
    $i=0;
    $ia=1;
    $ir=1;
    while(substr($line[$i],16,4) =~ /[BCDE234][A-Z][A-Z][A-Z]/) {
	$i++;
    }
    $atomName= substr($line[$i], 12,4); $atomName=~ s/^\s+|\s+$//g;
    if($atomName ne "CA"){next;}

    $resName1 = substr($line[$i], 17,3); $resName1 =~ s/^\s+|\s+$//g;
    $resNum1  = substr($line[$i], 22,4); $resNum1  =~ s/^\s+|\s+$//g;
    $chainID1 = substr($line[$i], 21,1); $chainID1 = "*" if($chainID1 eq " ");
    $x1       = substr($line[$i], 30,8); $x1       =~ s/^\s+|\s+$//g;
    $y1       = substr($line[$i], 38,8); $y1       =~ s/^\s+|\s+$//g;
    $z1       = substr($line[$i], 46,8); $z1       =~ s/^\s+|\s+$//g;
    printf OUT "ATOM  %5d  CA  %s %s%4d    %8.3f%8.3f%8.3f\n",$ia++,$resName1,$chainID1,$ir++,$x1,$y1,$z1;
    $i++;
    while(substr($line[$i],16,4) =~ /[BCDE234][A-Z][A-Z][A-Z]/) {
	$i++;
    }
    $resName = substr($line[$i], 17,3); $resName =~ s/^\s+|\s+$//g;
    $resNum  = substr($line[$i], 22,4); $resNum  =~ s/^\s+|\s+$//g;
    $chainID = substr($line[$i], 21,1); $chainID = "*" if($chainID eq " ");
    $x       = substr($line[$i], 30,8); $x       =~ s/^\s+|\s+$//g;
    $y       = substr($line[$i], 38,8); $y       =~ s/^\s+|\s+$//g;
    $z       = substr($line[$i], 46,8); $z       =~ s/^\s+|\s+$//g;
    printf OUT "ATOM  %5d  CA  %s %s%4d    %8.3f%8.3f%8.3f\n",$ia++,$resName,$chainID,$ir++,$x,$y,$z;
    $i++;
    for(;$i<$lineCnt;++$i) {
	if(substr($line[$i],16,4) =~ /[BCDE234][A-Z][A-Z][A-Z]/) {next;}

	$resName2 = substr($line[$i], 17,3); $resName2 =~ s/^\s+|\s+$//g;
	$resNum2  = substr($line[$i], 22,4); $resNum2  =~ s/^\s+|\s+$//g;
	$chainID2 = substr($line[$i], 21,1); $chainID2 = "*" if($chainID2 eq " ");
	$x2  = substr($line[$i], 30,8); $x2  =~ s/^\s+|\s+$//g;
	$y2  = substr($line[$i], 38,8); $y2  =~ s/^\s+|\s+$//g;
	$z2  = substr($line[$i], 46,8); $z2  =~ s/^\s+|\s+$//g;
# calculate C_b from Ca
	$xX=$x-$x1; $xY=$y-$y1; $xZ=$z-$z1;
	$yX=$x-$x2; $yY=$y-$y2; $yZ=$z-$z2;
	$zX=$xX+$yX; $zY=$xY+$yY; $zZ=$xZ+$yZ;
	$zP=sqrt($zX*$zX+$zY*$zY+$zZ*$zZ);
	$zX/=$zP; $zY/=$zP; $zZ/=$zP;
	$uX=$xY*$yZ-$xZ*$yY; $uY=$xZ*$yX-$xX*$yZ; $uZ=$xX*$yY-$xY*$yX;
	$uP=sqrt($uX*$uX+$uY*$uY+$uZ*$uZ);
	$uX/=$uP; $uY/=$uP; $uZ/=$uP;
	$vX=$zX * $R_sc[$Res{$resName}] * cos($pi * 37.6/180);
	$vY=$zY * $R_sc[$Res{$resName}] * cos($pi * 37.6/180);
	$vZ=$zZ * $R_sc[$Res{$resName}] * cos($pi * 37.6/180);
	$wX=$uX * $R_sc[$Res{$resName}] * sin($pi * 37.6/180);
	$wY=$uY * $R_sc[$Res{$resName}] * sin($pi * 37.6/180);
	$wZ=$uZ * $R_sc[$Res{$resName}] * sin($pi * 37.6/180);
	$CbX=$x+$vX+$wX; $CbY=$y+$vY+$wY; $CbZ=$z+$vZ+$wZ;
	if($resName ne "GLY"){
	    printf OUT "ATOM  %5d  SC  %s %s%4d    %8.3f%8.3f%8.3f\n",$ia++,$resName,$chainID,$ir-1,$CbX,$CbY,$CbZ;
	}
	printf OUT "ATOM  %5d  CA  %s %s%4d    %8.3f%8.3f%8.3f\n",$ia++,$resName2,$chainID2,$ir++,$x2,$y2,$z2;
	$resName1=$resName;	$resNum1=$resNum; $chainID1=$chainID;
	$x1=$x; $y1=$y; $z1=$z;
	$resName=$resName2;	$resNum=$resNum2; $chainID=$chainID2;
	$x=$x2; $y=$y2; $z=$z2;
    }
}

