#!/usr/bin/perl

# COUNTVIO
# Jane Caldwell 9/11/97
# Tabulates sizes and frequencies of constraint violations in X-PLOR.
# Accepts as input the output (log file) from any X-PLOR macro that 
# uses a "print" or "print threshold" command, e.g. refine.inp, 
# accept.inp, dgsa.inp.  Output from countvio is text in ascii format.
# This version supports the tabulation of violations of NOE, dihedral
# angle, and coupling constant restraints, as well as covalent bond and
# angle violations.

# Usage:  countvio.prl < xplor.log > violations.txt

# CAVEAT:
# Due to the algorithm used to search through the X-PLOR log file,
# extraneous comments on the same line as the "print" command, and
# comments containing the word "print" along with the words "angle",
# "bond", "noe", "improper", "cdih", or "coup" are to be avoided.

# Copyright (C) 1997, University of Wisconsin-Madison, all rights
# reserved.

# The program has been tested on an SGI Indigo running IRIX 5.3,
# and Perl version 4.0.
# Please send any comments to jane@nmrfam.wisc.edu.

# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.


$temp = <STDIN>;

$h = 0;
$i = 0;
$j = 0;
$k = 0;
$l = 0;
$m = 0;
$q = 0;
$strnum = 0;
$numviol = 0;
@noeres1[0] = 0;
@noeatm1[0] = 0;
@noeres2[0] = 0;
@noeatm2[0] = 0;
@noeviol[0] = 0;
@noefreq[0] = 0;
@dihres[0] = 0;
@dihatm[0] = 0;
@dihviol[0] = 0;
@dihfreq[0] = 0;
@coupres[0] = 0;
@coupatm[0] = 0;
@coupviol[0] = 0;
@coupfreq[0] = 0;
@bndres[0] = 0;
@bndatm[0] = 0;
@bndviol[0] = 0;
@bndfreq[0] = 0;
@angres[0] = 0;
@angatm[0] = 0;
@angviol[0] = 0;
@angfreq[0] = 0;
@impres[0] = 0;
@impatm[0] = 0;
@impviol[0] = 0;
@impfreq[0] = 0;
while ($temp ne "")
	{
	if ($temp =~ /prin/ && $temp =~ /noe/)
	   {
	   $strnum += 1;
		$temp = <STDIN>;
		$temp = <STDIN>;
		if ($temp =~ '\+\+\+\+\+')
		   {
		   &get_noe_viol;
		   $temp = <STDIN>;
		   while ($temp =~ '\=\=\=\=\=')
			{
			&get_noe_viol;
			$temp = <STDIN>;
			}
		   }
	   }
	if ($temp =~ /prin/ && $temp =~ /cdih/)
	   {
	   for ($count = 1; $count <=3; $count++)
		{
	   	$temp = <STDIN>;
		}
	   while ($temp =~ '\=\=\=\=\=')
		{
		&get_dihe_viol;
		$temp = <STDIN>
		}
	   }
	if ($temp =~ /prin/ && $temp =~ /coup/ && $temp =~ /thres/)
		{
		$temp = <STDIN>;
		$temp = <STDIN>;
		$temp = <STDIN>;
		$temp = <STDIN>;
		&get_coup_viol;
		}
	if ($temp =~ /prin/ && $temp =~ /bond/)
		{
		$temp = <STDIN>;
		$temp = <STDIN>;
		&get_bond_viol;
		}
	if ($temp =~ /prin/ && $temp =~ /angl/)
		{
		$temp = <STDIN>;
		$temp = <STDIN>;
		&get_angle_viol;
		}
	if ($temp =~ /prin/ && $temp =~ /impr/)
		{
		$temp = <STDIN>;
		$temp = <STDIN>;
		&get_improper_viol;
		}
	$temp = <STDIN>;
	}

$num_noe = $h;
$num_dihe = $j;
$num_coup = $q;
$num_bond = $k;
$num_angle = $l;
$num_improper = $m;

######################################
# PRINT OUT VIOLATIONS TO STDIN

if ($num_noe > 0)
	{
	&print_noe_viol;
	}
else 
	{ 
	print "\nNO NOE CONSTRAINT VIOLATIONS\n"; 
	}

if ($num_dihe > 0)
	{
	&print_dihe_viol;
	}
else 
	{ 
	print "\nNO DIHEDRAL ANGLE CONSTRAINT VIOLATIONS\n"; 
	}

if ($num_coup > 0)
	{
	&print_coup_viol;
	}
else 
	{ 
	print "\nNO J-COUPLING CONSTANT CONSTRAINT VIOLATIONS\n"; 
	}

if ($num_bond > 0)
	{
	&print_bond_viol;
	}
else 
	{ 
	print "\nNO BOND LENGTH VIOLATIONS\n"; 
	}

if ($num_angle > 0)
	{
	&print_angle_viol;
	}
else 
	{ 
	print "\nNO VALENCE ANGLE VIOLATIONS\n"; 
	}

if ($num_improper > 0)
	{
	&print_improper_viol;
	}
else 
	{ 
	print "\nNO IMPROPER ANGLE VIOLATIONS\n"; 
	}


#######################################
# SUBROUTINES

sub get_noe_viol
	{
                   {
                   until ($temp =~ 'set\-i\-atoms')
                        {
                        $temp = <STDIN>;
                        }
                   $temp = <STDIN>;
                   chop($temp);
                   @tmparr = split(/ +/,$temp);
	
                   $tmpres1 = @tmparr[2];
                   $tmpatm1 = @tmparr[4];
                   $temp = <STDIN>;
			
		   # special case: 1st part of constraint is pseudoatom
		   while ($temp !~ 'set\-j\-atoms')
                        {
                        @tmparr = split(/ +/,$temp);
			if ($tmpres1 != @tmparr[2])
				{ 
                        	$tmpres1 .= ", ".@tmparr[2];
				}
                        $tmpatm1 .= ", ".@tmparr[4];
			$temp = <STDIN>;	# advance to next line
                        }
			
                   if ($temp =~ 'set\-j\-atoms')
                        {
                        $temp = <STDIN>;
                        chop($temp);
                        @tmparr = split(/ +/,$temp);
                        $tmpres2 = @tmparr[2];
                        $tmpatm2 = @tmparr[4];
                        }
                   else
                        {
                        die "error in file format; 2nd part of NOE constraint not found";
                        }
                   $temp = <STDIN>;

		# special case:  2nd part of constraint is pseudoatom
		   while ($temp !~ Delta)
			{
			@tmparr = split(/ +/,$temp);
			if ($tmpres2 != @tmparr[2])
				{ 
                        	$tmpres2 .= ", ".@tmparr[2];
				}
			$tmpatm2 .= ", ".@tmparr[4];
			$temp = <STDIN>;	# advance to next line
			}

		# read in size of violation, etc.

                   if ($temp =~ Delta)
                        {
                        @delarr = split(/ +/,$temp);
                        $tmpviol = @delarr[9];
			$numviol +=1;
                        }
                   else
                        {
                        die "error in file format; Delta for NOE not found";
                        }

		# check to see if this constraint is violated multiple times
		   $i = 0;
		   $newitem = 1;
		# search for a match, stop when you find one
			do
			   {
			   if ($tmpres1 eq @noeres1[$i] && $tmpres2 eq @noeres2[$i] && $tmpatm1 eq @noeatm1[$i] && $tmpatm2 eq @noeatm2[$i])
				{
				$newitem = 0;
                       	 	@noefreq[$i] += 1;
                       	 	@noeviol[$i] += $tmpviol;
				}
			   $i += 1;
			   }
		   	until ($newitem eq "0" || $i > $h );

		# otherwise, for first occurrence of the constraint violation
		   if ($newitem > 0)
			{
                        @noeres1[$h] = $tmpres1;
                        @noeatm1[$h] = $tmpatm1;
                        @noeres2[$h] = $tmpres2;
                        @noeatm2[$h] = $tmpatm2;
			@noefreq[$h] = 1;
			@noeviol[$h] = $tmpviol;
			$h +=1;
			}
                   }
	}

sub get_dihe_viol
	{
	$tmpres = "";
	$tmpatm = "";
	for ($count=1; $count<=4; $count++)
	  {
	  $temp = <STDIN>;
	  chop($temp);
	  @tmparr = split(/ +/,$temp);
	  $tmpres .= @tmparr[2];
	  $tmpatm .= @tmparr[4];
	  if ($count < 4)
		{
		$tmpres .=",";
		$tmpatm .=",";
		}
	  }
	$temp = <STDIN>;
	chop($temp);
	@tmparr = split(/ +/,$temp);
	$tmpviol = @tmparr[10];

	# check to see if this constraint is violated multiple times
	$i = 0;
	$newitem = 1;
	# search for a match, stop when you find one
	do
		{
		if ($tmpres eq @dihres[$i] && $tmpatm eq @dihatm[$i])
			{
			$newitem = 0;
			@dihfreq[$i] += 1;
			@dihviol[$i] += $tmpviol;
			}
		$i += 1;
		}
	until ($newitem eq "0" || $i > $j );

	# otherwise, for first occurrence of the constraint violation
	if ($newitem > 0)
		{
		@dihres[$j] = $tmpres;
		@dihatm[$j] = $tmpatm;
		@dihfreq[$j] = 1;
		@dihviol[$j] = $tmpviol;
		$j +=1;
		}
	}

sub get_coup_viol
	{
	$temp = <STDIN>;

	# stop when you get to a line containing "RMS diff"
	# or when you reach the end of the file
	until ($temp =~ /RMS diff/ || $temp =~ /X-PLOR/)
		{
		chop($temp);
		@tmparr = split(/ +/,$temp);	
		$tmpres = @tmparr[2].",".@tmparr[6].",".@tmparr[10].",".@tmparr[14];
		$tmpatm = @tmparr[4].",".@tmparr[8].",".@tmparr[12].",".@tmparr[16];
		$tmpviol = @tmparr[19];

        	# check to see if this constraint is violated multiple times
        	$i = 0;
        	$newitem = 1;
        	# search for a match, stop when you find one
        	do
                	{
                	if ($tmpres eq @coupres[$i] && $tmpatm eq @coupatm[$i])
                        	{
                        	$newitem = 0;
                        	@coupfreq[$i] += 1;
                        	@coupviol[$i] += $tmpviol;
                        	}
                	$i += 1;
                	}
        	until ($newitem eq "0" || $i > $q );

        	# otherwise, for first occurrence of the constraint violation
        	if ($newitem > 0)
                	{
                	@coupres[$q] = $tmpres;
                	@coupatm[$q] = $tmpatm;
                	@coupfreq[$q] = 1;
                	@coupviol[$q] = $tmpviol;
                	$q +=1;
			}

		$temp = <STDIN>;
		}
	}

sub get_bond_viol
	{
        $temp = <STDIN>;

        # search for the pattern: ( xx...xx | xx...xx )
        while ($temp =~ /\(.+\|.+\)/)
                {
                chop($temp);
                @tmparr = split(/ +/,$temp);
                $tmpres = @tmparr[2].",".@tmparr[5];
                $tmpatm = @tmparr[3].",".@tmparr[6];
                $tmpviol = @tmparr[10];

                # check to see if this constraint is violated multiple times
                $i = 0;
                $newitem = 1;
                # search for a match, stop when you find one
                do
                        {
                        if ($tmpres eq @bndres[$i] && $tmpatm eq @bndatm[$i])
                                {
                                $newitem = 0;
                                @bndfreq[$i] += 1;
                                @bndviol[$i] += $tmpviol;
                                }
                        $i += 1;
                        }
                until ($newitem eq "0" || $i > $k );

                # otherwise, for first occurrence of the constraint violation
                if ($newitem > 0)
                        {
                        @bndres[$k] = $tmpres;
                        @bndatm[$k] = $tmpatm;
                        @bndfreq[$k] = 1;
                        @bndviol[$k] = $tmpviol;
                        $k +=1;
                        }

                $temp = <STDIN>;
                }
	}

sub get_angle_viol
	{
	$temp = <STDIN>;

	# search for the pattern: ( xx...xx | xx...xx | xx...xx )
	while ($temp =~ /\(.+\|.+\|.+\)/)
		{
		chop($temp);
		@tmparr = split(/ +/,$temp);	
		$tmpres = @tmparr[2].",".@tmparr[5].",".@tmparr[8];
		$tmpatm = @tmparr[3].",".@tmparr[6].",".@tmparr[9];
		$tmpviol = @tmparr[13];

        	# check to see if this constraint is violated multiple times
        	$i = 0;
        	$newitem = 1;
        	# search for a match, stop when you find one
        	do
                	{
                	if ($tmpres eq @angres[$i] && $tmpatm eq @angatm[$i])
                        	{
                        	$newitem = 0;
                        	@angfreq[$i] += 1;
                        	@angviol[$i] += $tmpviol;
                        	}
                	$i += 1;
                	}
        	until ($newitem eq "0" || $i > $l );

        	# otherwise, for first occurrence of the constraint violation
        	if ($newitem > 0)
                	{
                	@angres[$l] = $tmpres;
                	@angatm[$l] = $tmpatm;
                	@angfreq[$l] = 1;
                	@angviol[$l] = $tmpviol;
                	$l +=1;
			}

		$temp = <STDIN>;
		}
	}

sub get_improper_viol
	{
	$temp = <STDIN>;

	# search for the pattern: ( xx...xx | xx...xx | xx...xx )
	while ($temp =~ /\(.+\|.+\|.+\|.+\)/)
		{
		chop($temp);
		@tmparr = split(/ +/,$temp);	
		$tmpres = @tmparr[2].",".@tmparr[5].",".@tmparr[8].",".@tmparr[11];
		$tmpatm = @tmparr[3].",".@tmparr[6].",".@tmparr[9].",".@tmparr[12];
		$tmpviol = @tmparr[16];

        	# check to see if this constraint is violated multiple times
        	$i = 0;
        	$newitem = 1;
        	# search for a match, stop when you find one
        	do
                	{
                	if ($tmpres eq @impres[$i] && $tmpatm eq @impatm[$i])
                        	{
                        	$newitem = 0;
                        	@impfreq[$i] += 1;
                        	@impviol[$i] += $tmpviol;
                        	}
                	$i += 1;
                	}
        	until ($newitem eq "0" || $i > $m );

        	# otherwise, for first occurrence of the constraint violation
        	if ($newitem > 0)
                	{
                	@impres[$m] = $tmpres;
                	@impatm[$m] = $tmpatm;
                	@impfreq[$m] = 1;
                	@impviol[$m] = $tmpviol;
                	$m +=1;
			}

		$temp = <STDIN>;
		}
	}


sub print_noe_viol
	{
	$~ = "NOEHEADER";
	write;
	for ($count=0; $count < $num_noe; $count++)
		{
		&write_noe(@noeres1[$count],@noeatm1[$count],@noeres2[$count],
			@noeatm2[$count],@noeviol[$count],@noefreq[$count]);
		}
	}

sub print_dihe_viol
	{
	$~ = "DIHE_HEADER";
	write;
	for ($count=0; $count < $num_dihe; $count++)
		{
		&write_dihe(@dihres[$count],@dihatm[$count],@dihviol[$count],@dihfreq[$count]);
		}
	}

sub print_coup_viol
	{
	$~ = "COUP_HEADER";
	write;
	for ($count=0; $count < $num_coup; $count++)
		{
		&write_coup(@coupres[$count],@coupatm[$count],@coupviol[$count],@coupfreq[$count]);
		}
	}

sub print_bond_viol
	{
	$~ = "BONDHEADER";
	write;
	for ($count=0; $count < $num_bond; $count++)
		{
		&write_bond(@bndres[$count],@bndatm[$count],@bndviol[$count],@bndfreq[$count]);
		}
	}

sub print_angle_viol
	{
	$~ = "ANGLEHEADER";
	write;
	for ($count=0; $count < $num_angle; $count++)
		{
		&write_angle(@angres[$count],@angatm[$count],@angviol[$count],@angfreq[$count]);
		}
	}

sub print_improper_viol
	{
	$~ = "IMPROPERHEADER";
	write;
	for ($count=0; $count < $num_improper; $count++)
		{
		&write_improper(@impres[$count],@impatm[$count],@impviol[$count],@impfreq[$count]);
		}
	}

sub write_noe
	{
	local($res1,$atm1,$res2,$atm2,$viol,$freq) = @_;
	$~ = "WRITENOE";
	write;
	}

sub write_dihe
	{
	local($res,$atm,$viol,$freq) = @_;
	$~ = "WRITEDIHE";
	write;
	}

sub write_coup
	{
	local($res,$atm,$viol,$freq) = @_;
	$~ = "WRITECOUP";
	write;
	}

sub write_bond
	{
	local($res,$atm,$viol,$freq) = @_;
	$~ = "WRITEBOND";
	write;
	}

sub write_angle
	{
	local($res,$atm,$viol,$freq) = @_;
	$~ = "WRITEANGLE";
	write;
	}

sub write_improper
	{
	local($res,$atm,$viol,$freq) = @_;
	$~ = "WRITEIMPROPER";
	write;
	}

format NOEHEADER =

DISTANCE CONSTRAINT VIOLATIONS
Residue 1  Atom 1         Residue 2  Atom 2         Avg. Size (A)   No. Times   (Type)
======================================================================================
.

format WRITENOE =
@<<<<<     @<<<<<<<<<<<<< @<<<<<     @<<<<<<<<<<<<< @<<<<<          @<<<<       NOE
$res1,$atm1,$res2,$atm2,$viol/$freq,$freq
.

format DIHE_HEADER =

DIHEDRAL ANGLE CONSTRAINT VIOLATIONS
Residues          Atoms                             Avg. Size (deg) No. Times   (Type)
======================================================================================
.

format WRITEDIHE =
@<<<<<<<<<<<<     @<<<<<<<<<<<<<<<<<<<<<<<<<<<      @<<<<<          @<<<<       CDIH
$res,$atm,$viol/$freq,$freq
.

format COUP_HEADER =

J-COUPLING CONSTANT CONSTRAINT VIOLATIONS
Residues          Atoms                             Avg. Size (Hz)  No. Times   (Type)
======================================================================================
.

format WRITECOUP =
@<<<<<<<<<<<<     @<<<<<<<<<<<<<<<<<<<<<<<<<<<      @<<<<<          @<<<<       COUP
$res,$atm,$viol/$freq,$freq
.

format BONDHEADER =

BOND LENGTH VIOLATIONS
Residues          Atoms                             Avg. Size (A)   No. Times   (Type)
======================================================================================
.

format WRITEBOND =
@<<<<<<<<<<<<     @<<<<<<<<<<<<<<<<<<<<<<<<<<<      @<<<<<          @<<<<       BOND
$res,$atm,$viol/$freq,$freq
.

format ANGLEHEADER =

VALENCE ANGLE VIOLATIONS
Residues          Atoms                             Avg. Size (deg) No. Times   (Type)
======================================================================================
.

format WRITEANGLE =
@<<<<<<<<<<<<     @<<<<<<<<<<<<<<<<<<<<<<<<<<<      @<<<<<          @<<<<       ANGLE
$res,$atm,$viol/$freq,$freq
.

format IMPROPERHEADER =

IMPROPER ANGLE VIOLATIONS
Residues          Atoms                             Avg. Size (deg) No. Times   (Type)
======================================================================================
.

format WRITEIMPROPER =
@<<<<<<<<<<<<<<<< @<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< @<<<<<          @<<<<       IMPROP 
$res,$atm,$viol/$freq,$freq
.