#!/usr/bin/perl

#this script takes a *.pdb file as the first argument
#figures out volume and mw and psvol based on 
#method laid out in Harpaz et al, Structure, 15July94, pp641-649.
# editted from mk_char_seq_prot.pl

# 31Aug98  
# will now differentiate between sequence file and not

# form :: calc-psvol.pl [args val] *.(pdb/fasta)

# 8Sep98
# giving up on trying to get it to agree with previous paper
# will use pdb, to get correct disulfide data
# order:
# 1. finds disulfides from header
# 2. counts amino acids
# 3. figures mw and vol from
# 4. adjusts for disulfides and/or ends
# 5. outputs

# subroutine to return variable after command line argument
# usage: $stt = varafteroption("stt", 1.00, "start oxygen radii at ");
sub varafteroption {
  my $changed=0;
  for($ii=0;$ii<@ARGV;$ii++)  {
    if($ARGV[$ii] eq $_[0])  {
      $var = $ARGV[$ii+1]; $changed++;
      print STDERR "NEW $_[2]($_[0]) $var\n";
      last;
    }
  }
  if(!$changed)  {
    $var = $_[1];
    print STDERR "DEFAULT $_[2]($_[0]) $var\n";
  }
  return $var;
}

$volfile = varafteroption(vol,"CHC-res-vol", "volumes from file ");
open(VOL,"$volfile") || die "couldn't open $volfile\n$!\n";
while(<VOL>)  {
  next if /res/ || /#/; 
  chop; @tt=split' ';
  $vol{$tt[0]}=$tt[1];
  $svol{$tt[2]}=$tt[1];
}
close(VOL);

$mwfile = varafteroption(mw,"MW-res", "molecular weights from file ");
open(MW,"$mwfile") || die "couldn't open $mwfile\n$!\n";
while(<MW>)  {
  next if /res/ || /#/; 
  chop; @tt=split' ';
  $mw{$tt[0]}=$tt[1];
  $smw{$tt[2]}=$tt[1];
}
close(MW);

$doss = varafteroption(ss, 1, "doing ss ");

if (@ARGV>0)  { $nm=$ARGV[$#ARGV]; }
else { die "need single argument :: calc-psvol.pl [variables] filename\n"; }

#pdb file:
#ATOM      0  N   LYS A   1       3.240  10.040  10.380  1.00 -0.31
#ATOM      1  CA  LYS A   1       2.390  10.410   9.250  1.00 -0.20

$nres=0; $nlys=0; $narg=0; $nglu=0; $nasp=0;
$tmw=0; $tvol=0; $chains=0; $nss=0;
if ($nm =~ /pdb/)  {
  $pdbfile = $nm;
  open(PDB,$pdbfile) || die "couldn't open $pdbfile\n$!\n";

  while ($line=<PDB>)  {
    chop($line); @in=split(' ',$line);
    
    if($in[0] eq "SSBOND" && $doss)  {
      $nss++;
      if($in[3]=~/\D+/)  { $chains=1; } # don't really use
    }
    elsif($in[0] eq "ATOM" && $line =~ "CA")  {
      $res=$in[3]; $nres++;
      if($res eq "LYS") { $nlys++; }
      elsif($res eq "ARG") { $narg++; }
      elsif($res eq "GLU") { $nglu++; }
      elsif($res eq "ASP") { $nasp++; }
      $tmw+=$mw{$res}; $tvol+=$vol{$res};
    }
  }
  close(PDB);
}
elsif ($nm =~ /fasta/)  {
  open(OO,$nm) || die "couldn't open $pdbfile\n$!\n";

  while (<OO>)  {
    if (/^>/ .. /^$/)  {
      next if /^>/ || /^$/; chop; @in=split/ */;
      foreach $a (@in)  {
	$tmw+=$smw{$a}; $tvol+=$svol{$a};

	if($a eq "K") { $nlys++; }
	elsif($a eq "R") { $narg++; }
	elsif($a eq "D") { $nglu++; }
	elsif($a eq "E") { $nasp++; }

	$nres++;
      }
    }
  }
  close(OO);  
}

if(varafteroption("acetyl", 0, "including acetyl weights ")) { $tmw+=(24+16+3); }
$ex=varafteroption("ex", 0, "extra info ");

# calculate psvol
if($doss) {
    $dvds=$vol{"CYS"}-$vol{"CSS"};
    printf STDERR "CYS - CSS vol difference is %.2lf\n", $dvds;
    $tvol-=$nss*2*$dvds; 
    $dmds=$mw{"CYS"}-$mw{"CSS"};
    $tmw-=$nss*2*$dmds;
}				
# add extra oxygen at c-terminus and extra H at n-terminus
if(varafteroption("tt", 1, "including termini weights ")) { 
    $tmw+=(16+1); 
    if($vol{"Otm"}) {
	printf STDERR "using input %.2lf value for extra c terminal oxygen volume.\n", $vol{"Otm"};
	$tvol+=$vol{"Otm"};
    }
    else {
	printf STDERR "using default 18.0  value for extra c terminal oxygen volume.\n";
	$tvol+=18;
    }
    if($vol{"Htm"})  {
	printf STDERR "using input %.2lf value for extra n terminal hydrogen voluem.\n", $vol{"Htm"};
	$tvol+=$vol{"$Htm"};
    }
    else   {
	printf STDERR "not using value to correct of extra voluem at n terminus\n";
    }
}

# take 10 off for every carboxyl oxygen (two per asp and glu and two for end)
# take 18 off for amino and guanadino (include end)
if(varafteroption("el", 0, "taking off for electroconstriction ")) { 
    $tvol=$tvol-(2*10*($nasp+$nglu+1))-(18*($nlys+$narg+1));
}

printf "%15s res %d ss %d v %.2lf w %.2lf ", $nm, $nres, $nss, $tvol, $tmw;
if($ex) { printf "- D %d E %d _ K %d R %d ", $nasp, $nglu, $nlys, $narg; } 
$psvol=0.6023*$tvol/$tmw;
printf "psv %.3lf\n", $psvol;

__END__



