#!/usr/sbin/perl
# Add Na+ CIONS to a DNA PDB file
#
# usage: addcio.pl input.pdb outputname.pdb
#
#                                B.A. Luxon  04/05/98
$pdbfile = shift(@ARGV);
$outfile = shift(@ARGV);

open(PDB, "$pdbfile") || die "Cannot open file $pdbfile";
open(CIO, ">$outfile") || die "Cannot open file $outfile";

@pdb = <PDB>; # Read into array

$rmin = 99;
$rmax = 0;

foreach (@pdb) {
   next if (/Na/ || /CIO/ || /IB/ || /WAT/ || /IP/);
   if (/ATOM/) {
      chomp;
      print CIO "$_\n";
      
      if (/H3T/) {
         ($a,$anum,$anam,$rnam,$rnum,$x,$y,$z) = split('\s+',$_);      
          $amax = $anum;
          next;
      }

      if (/O.P/ || /S.P/ || / P /) {
         ($a,$anum,$anam,$rnam,$rnum,$x,$y,$z) = split('\s+',$_);
         /.1P/ && (($X1P[$rnum] = $x) && ($Y1P[$rnum] = $y) && ($Z1P[$rnum] = $z));
         /.2P/ && (($X2P[$rnum] = $x) && ($Y2P[$rnum] = $y) && ($Z2P[$rnum] = $z));
         / P / && (($XP[$rnum] = $x) && ($YP[$rnum] = $y) && ($ZP[$rnum] = $z));
         ($rnum < $rmin) && ($rmin = $rnum);
         ($rnum > $rmax) && ($rmax = $rnum);
         / P / && (push @Rnums,$rnum);
      }      
   }
}

foreach $rnum (@Rnums) {
   $rmax++ && $amax++;
   $Xc = $X1P[$rnum] + $X2P[$rnum] - $XP[$rnum];
   $Yc = $Y1P[$rnum] + $Y2P[$rnum] - $YP[$rnum];
   $Zc = $Z1P[$rnum] + $Z2P[$rnum] - $ZP[$rnum];
   write CIO;
}

#
format CIO =
ATOM  @#### Na+  Na+  @###    @###.###@###.###@###.### 
      $amax,          $rmax,  $Xc,    $Yc,    $Zc
.