#!/usr/local/bin/perl
# cdih_measure Last modified: 12-22-95  JPL
# Used to output tortion angles (heavy atom) from RNA/DNA .PDB files.
# Input file must be of "x-plor" std form:
# ATOM atom# atom_type res_type res# x y z

# Modified 22 Sep 96 by Dave Schweisguth <dcs@proton.chem.yale.edu>
# to parse PDB file based on columns, not whitespace

%types = (
	"alpha", ["O3'","P","O5'","C5'"],
	"beta", ["P","O5'","C5'","C4'"],
	"gamma",["O5'","C5'","C4'","C3'"],
	"chi", ["O4'","C1'","xx","xx"],
	"eps", ["C4'","C3'","O3'","P"],
	"zeta", ["C3'","O3'","P","O5'"],
	"nu0", ["C4'","O4'","C1'","C2'"],
	"nu1", ["O4'","C1'","C2'","C3'"],
	"nu2", ["C1'","C2'","C3'","C4'"],
	"nu3", ["C2'","C3'","C4'","O4'"],
	"nu4", ["C3'","C4'","O4'","C1'"]
);

$seg_num=0;

if ($ARGV[0] eq "define") {
	print "\n    Name = Definition\n";
        print "   alpha = O3'(n-1)-P-O5'-C5'            nu0 = C4'-O4'-C1'-C2'\n";
        print "    beta = P-O5'-C5'-C4'                 nu1 = O4'-C1'-C2'-C3'\n";
        print "   gamma = O5'-C5'-C4'-C3'               nu2 = C1'-C2'-C3'-C4'\n";
        print " epsilon = C4'-C3'-O3'-P(n+1)            nu3 = C2'-C3'-C4'-O4'\n";
        print "    zeta = C3'-O3'-P(n+1)-O5'(n+1)       nu4 = C3'-C4'-O4'-C1'\n";
        print "chi(pur) = O4'-C1'-N9-C4\n";
        print "chi(pyr) = O4'-C1'-N1-C2\n";
	print "\nStandard Values:";
	print "
segid  num res  alpha    beta   gamma     eps    zeta     chi     nu0     nu1     nu2     nu3     nu4
------ --- ---  ------ ------- ------- ------- ------- ------- ------- ------- ------- ------- -------
A-RNA           -62    -180      47    -152     -74    -166       3     -25      37     -36      21
B-DNA           -46.8  -146      36.4   155     -95.2   -97.8    -4.2    24.9   -34.9    33.3   -18.3
A-DNA           -83.9  -152.1    45.5    84.3   179.5   154.2     3.52  -26      37.1   -36      20.5\n\n";
	shift (@ARGV);
	}

foreach (<>) {

    # Can't do the match and assignment in one step, because we depend on the
    #   last assigned value of $num to be correct ... pity.

    next unless /^ATOM  .{5} (.{4}).(.{3}) .(.{4}).   (.{8})(.{8})(.{8}).{6}.{6} .{3}  (.{4})/;
    ($atom, $res, $num, $x, $y, $z, $segid) = ($1, $2, $3, $4, $5, $6, $7);

    # Trim whitespace from each field

    foreach $i ($atom, $res, $num, $x, $y, $z, $segid) {
	$i =~ s/^\s*([^\s]*)\s*$/$1/;
    }

    # Die if a field is empty which isn't allowed to be empty

    foreach $i ($atom, $res, $num, $x, $y, $z) {
	die "cdih_measure: Bad PDB record on input line ", __LINE__, ":\n$_"
	    if $i eq '';
    }
    
    if ($segid eq "") {
	$res[$num] = $res;
	$xyz[$num]{$atom} = [ $x, $y, $z ];
    } else {
	$segs=1;  # flag saying that there are named segments in the .pdb file.
	if ($segid ne $segid[$seg_num]) {
	    ++$seg_num;
	    $segid[$seg_num]=$segid;
	}
        $res[$num]{$segid} = $res;
        $xyz[$num]{$atom,$segid} = [ $x, $y, $z ];
    }
}

if ($seg_num eq 0) {$seg_num = 1;}

print "Number of segments = $seg_num\n";
print "Number of residues = $num\n";
print "segid  num res  alpha    beta   gamma     eps    zeta     chi     nu0     nu1     nu2     nu3     nu4\n";
print "------ --- ---  ------ ------- ------- ------- ------- ------- ------- ------- ------- ------- -------\n";

foreach $h (1 .. $seg_num) {
	$segid = $segid[$h];

foreach $i (1 .. $num) {
    foreach $j (sort keys %types) {
	if ($j eq "chi") {
	    ($atom_type1, $res_type1, $res_num1, $xyz1) =
		($types{$j}[0], $res[$i]{$segid}, $i, $xyz[$i]{$types{$j}[0],$segid});
	    ($atom_type2, $res_type2, $res_num2, $xyz2) =
		($types{$j}[1], $res[$i]{$segid}, $i, $xyz[$i]{$types{$j}[1],$segid});
	    if (($res[$i]{$segid} eq 'A') || ($res[$i]{$segid} eq 'G')) {
		($atom_type3, $res_type3, $res_num3, $xyz3) =
		    ("N9", $res[$i]{$segid}, $i, $xyz[$i]{"N9",$segid});
		($atom_type4, $res_type4, $res_num4, $xyz4) =
		    ("C4", $res[$i]{$segid}, $i, $xyz[$i]{"C4",$segid});
	    } else {
		($atom_type3, $res_type3, $res_num3, $xyz3) =
		    ("N1", $res[$i]{$segid}, $i, $xyz[$i]{"N1",$segid});
		($atom_type4, $res_type4, $res_num4, $xyz4) =
		    ("C2", $res[$i]{$segid}, $i, $xyz[$i]{"C2",$segid});
	    }
	} else {
	    if ($j eq "alpha") {
		if ($i == 1) {
		    $angle{'alpha'} = 0; next;
		} else {
		    ($atom_type1, $res_type1, $res_num1, $xyz1) =
		        ($types{$j}[0], $res[$i]{$segid}, $i-1, $xyz[$i-1]{$types{$j}[0],$segid});
		}
	    } else {
		($atom_type1, $res_type1, $res_num1, $xyz1) =
		    ($types{$j}[0], $res[$i]{$segid}, $i, $xyz[$i]{$types{$j}[0],$segid});
	    }

	    ($atom_type2, $res_type2, $res_num2, $xyz2) =
		($types{$j}[1], $res[$i]{$segid}, $i, $xyz[$i]{$types{$j}[1],$segid});

	    if ($j eq "zeta") {
		if ($res[$i+1]{$segid}) {
		    ($atom_type3, $res_type3, $res_num3, $xyz3) =
			($types{$j}[2], $res[$i+1]{$segid}, $i+1, $xyz[$i+1]{$types{$j}[2],$segid});
		} else {
		    $angle{'zeta'} = 0;
		    next;
		}
	    } else {
		($atom_type3, $res_type3, $res_num3, $xyz3) =
		    ($types{$j}[2], $res[$i], $i, $xyz[$i]{$types{$j}[2],$segid});
	    }

	    if (($j eq "zeta") && $res[$i+1]{$segid}) {
		($atom_type4, $res_type4, $res_num4, $xyz4) =
		    ($types{$j}[3], $res[$i+1]{$segid}, $i+1, $xyz[$i+1]{$types{$j}[3],$segid});
	    } elsif ($j eq "eps") {
		if ($res[$i+1]{$segid}) {
		     ($atom_type4, $res_type4, $res_num4, $xyz4) =
			($types{$j}[3], $res[$i+1]{$segid}, $i+1, $xyz[$i+1]{$types{$j}[3],$segid});
		} else {
		    $angle{'eps'} = 0;
		    next;
		}
	    } else {
		($atom_type4, $res_type4, $res_num4, $xyz4) =
		    ($types{$j}[3], $res[$i]{$segid}, $i, $xyz[$i]{$types{$j}[3],$segid});
	    }
	}

	($x1, $y1, $z1) = @$xyz1;
	($x2, $y2, $z2) = @$xyz2;
	($x3, $y3, $z3) = @$xyz3;
	($x4, $y4, $z4) = @$xyz4;

	$r21x=$x2-$x1; $r21y=$y2-$y1; $r21z=$z2-$z1;
	$r23x=$x2-$x3; $r23y=$y2-$y3; $r23z=$z2-$z3;
	$r34x=$x3-$x4; $r34y=$y3-$y4; $r34z=$z3-$z4;
	$Ax=($r21y*$r23z)-($r21z*$r23y);
	$Ay=($r21z*$r23x)-($r21x*$r23z);
	$Az=($r21x*$r23y)-($r21y*$r23x);
	$Bx=($r34y*$r23z)-($r34z*$r23y);
	$By=($r34z*$r23x)-($r34x*$r23z);
	$Bz=($r34x*$r23y)-($r34y*$r23x);
	$magA=sqrt($Ax**2+$Ay**2+$Az**2);
	$magB=sqrt($Bx**2+$By**2+$Bz**2);
	$dot=($Ax*$Bx+$Ay*$By+$Az*$Bz);
	$arccos=$dot/($magA*$magB);
	$bot=$arccos;
	$top=sqrt(1-$arccos**2);
	$theta=(180/3.1415927)*atan2($top,$bot);

	# new shit to calculate sign
	# calc theta2 between A and r34 because if <90, theta
	# should be positive, if >90 theta should be negative
	$mag34=sqrt($r34x**2+$r34y**2+$r34z**2);
	$dot2=($r34x*$Ax+$r34y*$Ay+$r34z*$Az);
	$arccos2=$dot2/($magA*$mag34);
	$theta2=(180/3.1415927)*atan2(sqrt(1-$arccos2**2),$arccos2);
	if ($theta2 > 90) {$theta=$theta*-1;}

	# Set the output variables
	$angle{$j} = $theta;
    }
    write;
}
print "\n";
}

################################################

format STDOUT =
@>>>>> @>> @>> @###.## @###.## @###.## @###.## @###.## @###.## @###.## @###.## @###.## @###.## @###.## 
$segid,$res_num2,$res_type2,$angle{'alpha'},$angle{'beta'},$angle{'gamma'},$angle{'eps'},$angle{'zeta'},$angle{'chi'},$angle{'nu0'},$angle{'nu1'},$angle{'nu2'},$angle{'nu3'},$angle{'nu4'}
.