#
# Module to calc energy and entropy contributions for mm_pbsa
#
# Holger Gohlke
#   18.10.2001
#
########################################################################

package mm_pbsa_calceneent;
require Exporter;

@ISA         = qw(Exporter);
@EXPORT      = qw(calc_energy_entropy);
@EXPORT_OK   = qw();
$VERSION     = 1.00;

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

use strict;

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

sub calc_energy_entropy(){
##########################

  # Parameters: \%GEN,\%DEC,\%DEL,\%MOL,\%PRO
  my $r_gen = shift;
  my $r_dec = shift;
  my $r_del = shift;
  my $r_mol = shift;
  my $r_pro = shift;

  print "\n";
  print "=>> Calculating energy / entropy contributions\n";

  my $suf_list = "";
  $suf_list .= "_com " if($r_gen->{"COMPLEX"});
  $suf_list .= "_rec " if($r_gen->{"RECEPTOR"});
  $suf_list .= "_lig " if($r_gen->{"LIGAND"});

  # Loop over COM/REC/LIG
  #######################
  my $suf;
  foreach $suf (split/ +/,$suf_list){
    # Generic names
    my $ncrd  = $r_gen->{"PATH"} . $r_gen->{"PREFIX"} . $suf . ".crd.";
    my $nbase = $r_gen->{"PREFIX"} . $suf;
    my $npdb  = $nbase . ".pdb.";
    my $npqr  = $nbase . ".pqr.";
    my $nlog  = $nbase . ".delphilog.";
    my $nxyz  = $nbase . ".xyz.";
    my $nmol  = $nbase . ".mslog.";
    my $nout  = $nbase . ".all.out";

    # Result file
    open(OUT,">$nout") || die("$nout not opened\n");
    print OUT "MM\n" if($r_gen->{"MM"});
    print OUT "GB\n" if($r_gen->{"GB"});
    print OUT "PB\n" if($r_gen->{"PB"});
    print OUT "MS\n" if($r_gen->{"GB"} || $r_gen->{"MS"}); # Impicit use of sander SAS calc 
    print OUT "NM\n" if($r_gen->{"NM"});                   #   for GB if not MS set

    # Parmtop file
    my $npar;
    $npar = $r_gen->{"COMPT"} if($suf eq "_com");
    $npar = $r_gen->{"RECPT"} if($suf eq "_rec");
    $npar = $r_gen->{"LIGPT"} if($suf eq "_lig");

    # Generic input files
    my $nsan  = "sander" . $suf . ".in";
    my $nscl  = "sanmin" . $suf . ".in";
    my $nnmo  = "nmode" . $suf . ".in";

    if($r_gen->{"NM"}){
      if(! -e $nscl){
        die("    $nscl not found\n");
      }
      if(! -e $nnmo){
        die("    $nnmo not found\n");
      }
    }
    if(! -e $npar){
      die("$npar not found\n");
    }

    # Loop over all coordinates
    ###########################
    my $number = 1;
    while(1){
      my $ncrdfile = $ncrd . "$number";
      last if(! -e $ncrdfile);
      print "    Calc contrib for $ncrdfile\n";
      print OUT $number,"\n";

      # Calc MM and GB and SAS (if implicit)
      ######################################
      if($r_gen->{"MM"} || $r_gen->{"GB"}){
        my $sanout = "sander" . $suf . "." . $number . ".out";
        &calc_MM_GB_SAS(*OUT,$nsan,$sanout,$ncrdfile,$npar,$r_del,$r_gen,$r_pro);
      }

      # Calc entropy
      ##############
      if($r_gen->{"NM"}){
        my $sanout = "sanmin" . $suf . "." . $number . ".out";
        my $sanres = "sanmin" . $suf . "." . $number . ".restrt";
        my $nmodeout = "nmode" . $suf . "." . $number . ".out";
        &calc_NM(*OUT,$nscl,$nnmo,$ncrdfile,$npar,$sanout,$sanres,$nmodeout,$r_pro);
      }

      # Generate pdb file for PB or MS calc
      #####################################
      my $pdb = "";
      if(($r_gen->{"PB"} > 0 && $r_del->{"PROC"} == 1) || $r_gen->{"MS"} > 0){
        $pdb = $npdb . $number;
        &generate_pdb($pdb,$ncrdfile,$npar,$r_pro);
        &center_pdb($pdb);
      }

      # Calc PB
      #########
      if($r_gen->{"PB"}){
	if($r_del->{"PROC"} == 1){
	  my $pdbbase = $npdb . $number;
	  my $outbase = $nlog . $number;
	  &calc_delphi(*OUT,$pdbbase,$outbase,$r_del,$r_pro);
	}
	elsif($r_del->{"PROC"} == 2){
	  my $pbsaout = "pbsa" . $suf . "." . $number . ".out";
	  &calc_pbsa(*OUT,"pbsa.in",$pbsaout,$ncrdfile,$npar,$r_pro);
	}
      }

      # Calc MS
      ###########
      if($r_gen->{"MS"}){
        my $pqr = $npqr . $number;
        &generate_pqr($pdb,$pqr,$r_mol);

        my $mol = $nmol . $number;
        &calc_ms(*OUT,$pqr,$mol,$r_mol,$r_pro);
      }

      # Final clean up's
      if(-e "${npdb}${number}"){
        unlink $npdb . $number;
      }

      # Increase counter
      $number++;

    } # end while
    close(OUT);
  } # end for
}

sub calc_MM_GB_SAS(){
#####################

  # Parameters: *OUT,$nsan,$sanout,$ncrdfile,$npar,$r_del,$r_gen,$r_pro
  my $OUT = shift;
  my $nsan = shift;
  my $sanout = shift;
  my $ncrdfile = shift;
  my $npar = shift;
  my $r_del = shift;
  my $r_gen = shift;
  my $r_pro = shift;

  print "        Calc MM/GB/SAS\n";

  my $command = $r_pro->{"SANDER"} . " -O"   .
                                     " -i "  . $nsan .
                                     " -o "  . $sanout .
                                     " -c "  . $ncrdfile .
                                     " -p "  . $npar;

  system($command) && die("        $command not successful\n");
  
  # Write MM results
  open(IN,"$sanout") || die("        $sanout not opened\n");
  my $finalflg = 0;
  my $line;
  while(defined($line = <IN>)){
    if($r_gen->{"DC"}){
      if($line =~ /^TDC/ || $line =~ /^SDC/ || $line =~ /^BDC/){
        print OUT $line;
      }
    }
    else{
      if($line =~ /FINAL RESULTS/){
        $finalflg = 1;
      }
      elsif($finalflg > 0 &&
            ($line =~ /BOND +=/    ||
             $line =~ /ANGLE +=/   ||
             $line =~ /DIHED +=/   ||
             $line =~ /1-4 VDW +=/ ||
             $line =~ /1-4 EEL +=/ ||
             $line =~ /VDWAALS +=/ ||
             $line =~ /EEL +=/     ||
             ($line =~ /EGB +=/ && $r_gen->{"GB"} > 0) ||
             ($line =~ /ESURF +=/ && $r_gen->{"GB"} > 0 && $r_gen->{"MS"} == 0))
           ){
        $line =~ s/ ESURF +=/surface area =/; # Equal output as for MS
        print OUT $line;
      }
    }
  }
  close(IN);

  # Clean up MM
  unlink $sanout;
}

sub calc_pbsa(){
################

  # Parameters: $OUT, $pbsain, $pbsaout, $ncrdfile, $npar, $r_pro

  my $OUT = shift;
  my $pbsain = shift;
  my $pbsaout = shift;
  my $ncrdfile = shift;
  my $npar = shift;
  my $r_pro = shift;

  print "        Calc PBSA\n";

  my $command = $r_pro->{"PBSA"} . " -O"   .
                                   " -i "  . $pbsain .
                                   " -o "  . $pbsaout .
                                   " -c "  . $ncrdfile .
                                   " -p "  . $npar;

  system($command) && die("        $command not successful\n");
  
  # Write MM results
  open(IN,"$pbsaout") || die("        $pbsaout not opened\n");
  my $finalflg = 0;
  my $line;
  while(defined($line = <IN>)){
    if($line =~ /FINAL RESULTS/){
      $finalflg = 1;
    }
    elsif($finalflg > 0 && $line =~ /VDWAALS.+EEL += +(-?\d+\.\d+).+ERF/){
      # Convert to kT, because back-conversion to kcal/mol is done in 
      #   mm_pbsa_statistics, subroutine treat_special (to be compatible with delphi calculation)
      my $kcal2kt = 1000.0 * 4.184 / (300.0 * 8.314);
      my $rfe = $1 * $kcal2kt;
      printf OUT "%s %15.6f\n", "corrected reaction field energy:", $rfe;
    }
  }
  close(IN);

  # Clean up MM
  unlink $pbsaout;
}

sub calc_delphi(){
##################

  # Parameters: $OUT,$pdbbase,$outbase,\%DEL,\%PRO
  my $OUT = shift;
  my $pdbbase = shift;
  my $outbase = shift;
  my $r_del = shift;
  my $r_pro = shift;

  print "        Calc delphi\n";

  # Delphi calculation

  # Input PDB file
  unlink "fort.13" if(-l "fort.13");
  unlink "fort.13" if(-e "fort.13");
  symlink $pdbbase, "fort.13";

  # Do delphi (with focussing)
  my $multruns = 1;
  if((exists($r_del->{"SALT"}) && ($r_del->{"SALT"} > 0.0)) ||
     $r_del->{"REFE"} > 0){
    $multruns = 2;
  }
  my $focruns = $r_del->{"FOCUS"} + 1;

  my $j;
  for($j = 0; $j < $multruns; $j++){
    my $i;
    my $delphiout = "";
    for($i = 0; $i < $focruns; $i++){
      print "            Run $i $j\n";
      unlink "ARCDAT" if(-e "ARCDAT");
      unlink $delphiout if(-e $delphiout); # Remove delphiout from last run
      $delphiout = $outbase . "." . $i . "." . $j;
      my $delphiprm = "delphi.prm." . $i . "." . $j;
      my $command = $r_pro->{"DELPHI"} . " " . $delphiprm . " > " . $delphiout;
      system($command) && die("        $command not running properly\n");
    }
  }

  # Write PB results
  my $solv = 0.0;
  for($j = 0; $j < $multruns; $j++){
    my $delphiout = $outbase . "." . $r_del->{"FOCUS"} . "." . $j;
    open(IN,"$delphiout") || die("        $delphiout not opened\n");
    my $line;
    while(defined($line = <IN>)){
      if($r_del->{"REFE"} == 0){
	if($j == 0 && $line =~ /corrected reaction field energy: +(-?\d+\.\d+)/){
	  $solv += $1;
	}
	if($multruns == 2 && $line =~ /total energy \(including grid energy\): +(-?\d+\.\d+)/){
	  # solv = solv(SALT = 0.0) + total(SALT > 0.0) - total(SALT = 0.0)
	  #        delphi run 1       delphi run 2        delphi run 1
	  if($j == 0){
	    $solv -= $1;
	  }
	  else{
	    $solv += $1;
	  }
	}
      }
      else{ 
	# REFE > 0 does not work with SALT > 0.0 yet
	if($j == 0 && $line =~ /corrected reaction field energy: +(-?\d+\.\d+)/){
	  $solv -= $1; # 1.0/INDI calc
	}
	elsif($line =~ /corrected reaction field energy: +(-?\d+\.\d+)/){
	  $solv += $1; # EXDI/INDI calc
	}
      }
    }
    close(IN);
    unlink $delphiout;
  }
  printf OUT "%s %15.6f\n", "corrected reaction field energy:", $solv; 

  # Clean up PB
  unlink "fort.13";
  unlink "fort.14";
  unlink "ARCDAT";

  return;
}

sub calc_NM(){
##############

  # Parameters: *OUT,$nscl,$nnmo,$ncrdfile,$npar,$sanout,$sanres,$nmodeout,\%PRO
  my $OUT = shift;
  my $nscl = shift;
  my $nnmo = shift;
  my $ncrdfile = shift;
  my $npar = shift;
  my $sanout = shift;
  my $sanres = shift;
  my $nmodeout = shift;
  my $r_pro = shift;

  print "        Minimize structure\n";
  my $command = $r_pro->{"SANDER"} . " -O"  .
                                     " -i " . $nscl .
                                     " -o " . $sanout .
                                     " -c " . $ncrdfile .
                                     " -p " . $npar .
                                     " -r " . $sanres; 
  system($command) && die("        $command not running properly\n");

  print "        Calc entropy\n";
  $command = $r_pro->{"NMODE"} . " -O"  .
                                 " -i " . $nnmo .
                                 " -o " . $nmodeout .
                                 " -c " . $sanres .
                                 " -p " . $npar;
  system($command) && die("        $command not running properly\n");
       
  open(IN,"$nmodeout") || die("         $nmodeout not opened\n");
  my $line;
  while(defined($line = <IN>)){
    ###if($line =~ /(Total|translational|rotational|vibrational) +\d+\.\d+ +\d+\.\d+ +(\d+\.\d+)/){
    if($line =~ /(Total|translational|rotational|vibrational) +\S+ +\S+ +(\d+\.\d+)/){
      print OUT "$1 $2\n";
    }
  }
  close(IN);

  # Clean up
  unlink $sanout;
  unlink $sanres;
  unlink $nmodeout;
}

sub calc_ms(){
##############

  # Parameters: $out,$pqr,$mol,$r_mol,$r_pro
  my $OUT = shift;
  my $pqr = shift;
  my $mol = shift;
  my $r_mol = shift;
  my $r_pro = shift;

  print "        Calc MS\n";
  my $command = $r_pro->{"MS"} . " $pqr " . $r_mol->{"PROBE"} . " > " . $mol;
  system($command);

  open(IN,"$mol") || die("         $mol not opened\n");
  my $line;
  while(defined($line = <IN>)){
    if($line =~ /surface area =/){
      print OUT $line;
    }
  }
  close(IN);

  unlink $pqr;
  unlink $mol;
}

sub generate_pqr(){
###################

  # This implementation uses bondi radii

  # Parameters: $pdb,$pqr,$r_mol
  my $pdb = shift;
  my $pqr = shift;
  my $r_mol = shift;

  # Bondi radii + 1.4A and probe radius of 0.0A yields SAS
  # Bondi radii + 0.0A and probe radius of 1.4A yields molecular surface
  # Bondi radii + 0.0A and probe radius of 0.0A yields vdW surface
  my %exp_rad = (
		 "N"    => 1.550 + 1.400,
		 "H"    => 1.200 + 1.400,
		 "C"    => 1.700 + 1.400,
		 "O"    => 1.500 + 1.400,
		 "F"    => 1.470 + 1.400,
                 "P"    => 1.800 + 1.400,
		 "S"    => 1.800 + 1.400,
		 "FE"   => 1.300 + 1.400,
		 "Na+"  => 1.200 + 1.400,
		 "Cl-"  => 1.700 + 1.400,
		 "MG"   => 1.180 + 1.400,
		);

  print "        Generate PQR\n";
  make_pqr_file($pdb,$pqr,\%exp_rad);
}

sub make_pqr_file(){
####################

  # Parameters: $pdb,$pqr,\%exp_rad
  my $pdb = shift;
  my $pqr = shift;
  my $r_exp_rad = shift;

  open(PDB,"$pdb")  || die("        $pdb not opened\n");
  open(PQR,">$pqr") || die("        $pqr not opened\n");
  my $line;
  while(defined($line = <PDB>)){
    chomp($line);
    if($line =~ /(ATOM|HETATM)/){
      my(        $card, $atnum, $atm1, $atm2, $alt, $resname, $resno, $x, $y, $z) =
	 unpack("a6     a5 x    a      a3     a     a3 x2     a4 x4   a8  a8  a8", $line);

      my $atm = $atm1 . $atm2;
      my $atmtmp = $atm;  # Remove leading and trailing white spaces because patterns below
      $atmtmp =~ s/\s//g; #   are of the type /^\S+$/
 
      my $rad;
      if(exists $r_exp_rad->{"$atmtmp"}){
	# First check full atom name
	$rad = $r_exp_rad->{"$atmtmp"};
      }
      else{
	# Then check first character only
	$atmtmp = substr $atmtmp, 0, 1;
	if(exists $r_exp_rad->{"$atmtmp"}){
	  $rad = $r_exp_rad->{"$atmtmp"};
	}
	else{
	  print "        No radius found for $atm $atnum in residue $resname $resno\n";
	  die();
	}
      }

      $line = sprintf("%-6s%5d %1s%-3s%1s%-3s  %4d    %8.3f%8.3f%8.3f%6.2f%6.2f",
		      $card,
		           $atnum,
		               $atm1, $atm2,
                                      $alt,
                                         $resname,
                                               $resno,$x,  $y,  $z,  "00.000", $rad);

    } # endif
    print PQR $line,"\n";
  } # endwhile
  close(PQR);
  close(PDB);
}

sub generate_pdb(){
###################

  # Parameters: $pdb,$ncrdfile,$npar,$r_pro
  my $pdb = shift;
  my $ncrdfile = shift;
  my $npar = shift;
  my $r_pro = shift;

  print "        Generate PDB\n";

  unlink "tmp.pdb" if(-e "tmp.pdb");
  my $command = $r_pro->{"AMBPDB"} . " -p " . $npar .
                                     " < " . $ncrdfile . " > tmp.pdb 2> /dev/null";
  system("$command") && die("        $command not successful\n");
  
  # ----- The following implements the former "check_names" program
  #       However, ambpdb produces "better" pdb files than anal
  #         -> this part should be redundant.
  open(IN,"tmp.pdb") || die("        tmp.pdb not opened\n");
  open(PDB,">$pdb")  || die("        $pdb not opened\n");
  my $line;
  while(defined($line = <IN>)){
    #              ATOM            1000       1   HB   rest_of_line
    #              $1              $2     $3  $4  $5   $6
    chomp($line);
    if($line =~ /(ATOM +|HETATM +)(\d+)( +)( |\d)(\S+)( +.+$)/){
       $line = $1 . $2 . $3 . $5 . $4 . $6;
    }
    print PDB $line,"\n";
  }
  close(PDB);
  close(IN);
  # -----

  unlink "tmp.pdb" if(-e "tmp.pdb");
}

sub center_pdb(){
#################

  # Parameters: $pdb
  my $pdb = shift;

  print "        Center PDB\n";

  unlink "tmp.pdb" if(-e "tmp.pdb");

  open(IN,"$pdb") || die("        $pdb not opened\n");
  my @cont = <IN>;
  close(IN);

  my $line;
  my ($x,$y,$z);
  my $xcent = 0.0;
  my $ycent = 0.0;
  my $zcent = 0.0;
  my $count = 0;
  foreach $line (@cont){
    chomp($line);
                #ATOM      1 N    MET     1      55.891  31.265  43.091
    if($line =~ /ATOM(.{26})( *-?\d+\.\d+)( *-?\d+\.\d+)( *-?\d+\.\d+)/){
      $x = substr $line, 30, 8;
      $y = substr $line, 38, 8;
      $z = substr $line, 46, 8;
      $xcent += $x;
      $ycent += $y;
      $zcent += $z;
      $count++;
    }
    elsif($line =~ /ATOM/){
      die("        Unrecognized ATOM format in pdb: $line\n");
    }
  }

  $xcent /= $count;
  $ycent /= $count;
  $zcent /= $count;

  open(TMP,">tmp.pdb") || die("        tmp.pdb not opened\n");
  foreach $line (@cont){
    chomp($line);
                #ATOM      1 N    MET     1      55.891  31.265  43.091
    if($line =~ /ATOM(.{26})( *-?\d+\.\d+)( *-?\d+\.\d+)( *-?\d+\.\d+)/){
      $x = substr $line, 30, 8;
      $y = substr $line, 38, 8;
      $z = substr $line, 46, 8;
      $x -= $xcent;
      $y -= $ycent;
      $z -= $zcent;
      printf TMP "ATOM%s%8.2f%8.2f%8.2f\n", $1, $x, $y, $z;
    }
    else{
      print TMP $line,"\n";
    }
  }
  close(TMP);

  rename "tmp.pdb", $pdb;
}

1; # Necessary for package function