#!/usr/bin/perl
#removes dummy atoms around the heavy atoms
#arp/warp should not remove heavy atoms from models in later versions and then we don't have to do this

$tresh = 2.;

#open(HEAVY,"heavy.pdb_exp") || die "Problem opening HEAVY file in remove_dummy.";
#open(ARPIN,"arp_0.pdb") || die "Problem opening ARPIN file in remove_dummy.";
#open(ARPOUT,">arp_0.pdb_") || die "Problem opening ARPOUT file in remove_dummy.";
open(HEAVY,$ARGV[0]) || die "Problem opening HEAVY file ($ARGV[0]) in remove_dummy.";
open(ARPIN,$ARGV[1]) || die "Problem opening ARPIN file ($ARGV[1]) in remove_dummy.";
open(ARPOUT,">$ARGV[1]_") || die "Problem opening ARPOUT file in remove_dummy.";

#get the transformation matrix and its inverse
while ( $line=<ARPIN> ) 
{
  if ( $line =~ /SCALE1\s+(\d*\.*\d*)\s+(-?\d*\.*\d*)\s+(-?\d*\.*\d*)/ ) { $mat1[1]=$1; $mat1[2]=$2; $mat1[3]=$3; }
  if ( $line =~ /SCALE2\s+(\d*\.*\d*)\s+(\d*\.*\d*)\s+(-?\d*\.*\d*)/ ) { $mat2[1]=$1; $mat2[2]=$2; $mat2[3]=$3; }
  if ( $line =~ /SCALE3\s+(\d*\.*\d*)\s+(\d*\.*\d*)\s+(\d*\.*\d*)/ ) { $mat3[1]=$1; $mat3[2]=$2; $mat3[3]=$3; }
}
#print "$mat1[1] $mat2[2] $mat3[3]\n";
if ( $mat1[1] == 0. || $mat2[2] == 0. || $mat3[3] == 0.) { die "Problem reading ARPIN file in remove_dummy.";}
$invmat3[3] = 1./$mat3[3];
$invmat2[2] = 1./$mat2[2]; 
$invmat2[3] =  -$mat2[3]*$invmat2[2]*$invmat3[3];
$invmat1[1] = 1./$mat1[1];
$invmat1[2] = -$mat1[2]*$invmat1[1]*$invmat2[2];
$invmat1[3] = -$mat1[2]*$invmat1[1]*$invmat2[3] - $mat1[3]*$invmat1[1]*$invmat3[3];

close(ARPIN);
#open(ARPIN,"arp_0.pdb") || die "Problem opening ARPIN file 2 in remove_dummy.";
open(ARPIN,$ARGV[1]) || die "Problem opening ARPIN file 2 in remove_dummy.";


#read heavy atoms
$num_heavy=0;
while ( $line=<HEAVY> ) 
{
  if ( $line =~ /(HETATM|ATOM).{22}\s*(-*\d+\.*\d*)\s*(-*\d+\.*\d*)\s*(-*\d+\.*\d*)\s*(-*\d+\.*\d*)/ ) { 
	if ( $5 > 0.15 ) {
	  $heavy1[$num_heavy]=$2; 	$heavy2[$num_heavy]=$3; 	$heavy3[$num_heavy]=$4;
	  $num_heavy++;
	}
  }
  if ( $line =~ /HEAVY\s*\S+\s*(-*\d+\.*\d*)\s*(-*\d+\.*\d*)\s*(-*\d+\.*\d*)\s*(-*\d+\.*\d*)/ ) {
	if ( $4 > 0.15 ) {
	  $heavy1[$num_heavy]=$1; 	$heavy2[$num_heavy]=$2; 	$heavy3[$num_heavy]=$3;
	  $num_heavy++;
	}
  }
}
#for ($i=0; $i<$num_heavy; $i++) { print "$heavy1[$i] $heavy2[$i] $heavy3[$i]\n"; }
#move heavy atoms to defined unit cell
for ($h=0; $h<$num_heavy; $h++) 
{ 
  $fx = $mat1[1]*$heavy1[$h] + $mat1[2]*$heavy2[$h] + $mat1[3]*$heavy3[$h];
  $fy = $mat2[2]*$heavy2[$h] + $mat2[3]*$heavy3[$h];
  $fz = $mat3[3]*$heavy3[$h];

  while ( $fx<0. ) { $fx += 1.; }
  while ( $fy<0. ) { $fy += 1.; }
  while ( $fz<0. ) { $fz += 1.; }
  while ( $fx>1. ) { $fx -= 1.; }
  while ( $fy>1. ) { $fy -= 1.; }
  while ( $fz>1. ) { $fz -= 1.; }

  $heavy1[$h] = $invmat1[1]*$fx + $invmat1[2]*$fy + $invmat1[3]*$fz;
  $heavy2[$h] = $invmat2[2]*$fy + $invmat2[3]*$fz;
  $heavy3[$h] = $invmat3[3]*$fz;
#  print "HEAVY $h $heavy1[$h] $heavy2[$h] $heavy3[$h]\n";
}

#read dummy atoms and remove if too close to heavy atom
$num_rem = 0;
while ( $line=<ARPIN> ) 
{
  $write=1;
  if ( $line =~ /DUM  DUM.{8}\s*(-*\d+\.*\d*)\s*(-*\d+\.*\d*)\s*(-*\d+\.*\d*)/ ) {
#	print "$1 $2 $3\n";
	$dum1 = $1;	$dum2 = $2;	$dum3 = $3;

	#transform to defined unit cell
	$fx = $mat1[1]*$dum1 + $mat1[2]*$dum2 + $mat1[3]*$dum3;
	$fy = $mat2[2]*$dum2 + $mat2[3]*$dum3;
	$fz = $mat3[3]*$dum3;
	while ( $fx<0. ) { $fx += 1.; }
	while ( $fy<0. ) { $fy += 1.; }
	while ( $fz<0. ) { $fz += 1.; }
	while ( $fx>1. ) { $fx -= 1.; }
	while ( $fy>1. ) { $fy -= 1.; }
	while ( $fz>1. ) { $fz -= 1.; }
	$dum1 = $invmat1[1]*$fx + $invmat1[2]*$fy + $invmat1[3]*$fz;
	$dum2 = $invmat2[2]*$fy + $invmat2[3]*$fz;
	$dum3 = $invmat3[3]*$fz;

	#check distance
	for ($h=0; $h<$num_heavy; $h++) {
	  if ( ($dist=sqrt( ($dum1-$heavy1[$h])**2 + ($dum2-$heavy2[$h])**2 + ($dum3-$heavy3[$h])**2 )) < $tresh ) 
	  { 
		$write = 0; 
#		 print "$h  $dist      $line"; 
		$num_rem++;
	  }
	}
  }
  if ($write) { syswrite(ARPOUT,$line); }
}

close(HEAVY);
close(ARPIN);
close(ARPOUT);

#system("cp -f $ARGV[1] ${ARGV[1]}__");
#system("cp -f ${ARGV[1]}_ $ARGV[1]");
system("mv -f ${ARGV[1]}_ $ARGV[1]");
#print "number: $num_heavy\n";
print "Removed $num_rem dummy atoms too close to heavy atoms.\n";
