#
# Module with functions to create input for mm_pbsa "helper" programs
#
# Holger Gohlke
#   17.04.2002
#
########################################################################

package mm_pbsa_createinput;
require Exporter;

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

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

use strict;
use lib ("$ENV{AMBERHOME}/src/mm_pbsa");
use mm_pbsa_global qw(@GC_PAR);

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

sub create_input(){
###################

  # Parameters: \%GEN,\%DEC,\%SAN,\%DEL,\%GBO,\%MOL,\%NMO,\%MAK,\%PRO,\$TRA
  my $r_gen = shift;
  my $r_dec = shift;
  my $r_san = shift;
  my $r_del = shift;
  my $r_gbo = shift;
  my $r_mol = shift;
  my $r_nmo = shift;
  my $r_mak = shift;
  my $r_pro = shift;
  my $r_tra = shift;

  print "\n";
  print "=>> Creating input\n";

  if($r_gen->{"MM"} || $r_gen->{"GB"}){
    &create_sander_input($r_gen,$r_dec,$r_san,$r_del,$r_gbo);
  }

  if($r_gen->{"PB"}){
    if($r_del->{"PROC"} == 1){
      &create_delphi_input($r_del);
    }
    elsif($r_del->{"PROC"} == 2){
      &create_pbsa_input($r_gen, $r_dec, $r_san, $r_del);
    }
    elsif($r_del->{"PROC"} == 3){
      &create_iapbs_input($r_gen, $r_san, $r_del);
    }
  }

  if($r_gen->{"GC"} || $r_gen->{"AS"}){
    &create_makecrd_input($r_gen,$r_mak,$r_tra);
  }

  if($r_gen->{"NM"}){
    &create_nmode_input($r_gen,$r_nmo);
  }
  
  return;
}

sub create_nmode_input(){
#########################

  # Parameters: \%GEN,\%NMO
  my $r_gen = shift;
  my $r_nmo = shift;

  print "    Nmode input\n";

  my ($dielc, $maxcyc, $drms, $nmdrms);
  $dielc = $r_nmo->{"DIELC"};
  $maxcyc = $r_nmo->{"MAXCYC"};
  $drms = $r_nmo->{"DRMS"};
  $nmdrms = $drms * 10.0; # Increased drms for nmode due to differences 
                          # in gradient calc. with respect to sander.

  my $suf_list = "";
  $suf_list .= "_com " if($r_gen->{"COMPLEX"});
  $suf_list .= "_rec " if($r_gen->{"RECEPTOR"});
  $suf_list .= "_lig " if($r_gen->{"LIGAND"});
  my $suf;
  foreach $suf (split/ +/,$suf_list){
    my $name = "sanmin" . $suf . ".in";
    open(OUT,">$name") || die("    $name not opened\n");
    print OUT "File generated by mm_pbsa.pl\n";
    print OUT " &cntrl\n";
    ### print OUT "  ntxo   = 0,\n"; # Binary format is no longer supported
    print OUT "  ntf    = 1,       ntb    = 0,\n";
    print OUT "  dielc  = $dielc,\n";
    print OUT "  cut    = 99.0,    nsnb   = 99999,\n";
    print OUT "  scnb   = 2.0,     scee   = 1.2,\n";
    print OUT "  imin   = 1,       maxcyc = $maxcyc,\n";
    print OUT "  ncyc   = 0,       drms   = $drms\n";
    print OUT " &end\n";
    print OUT " &ewald\n";
    print OUT "  eedmeth= 5,\n";
    print OUT " &end\n";
    close(OUT);

    $name = "nmode" . $suf . ".in";
    open(OUT,">$name") || die("    $name not opened\n");
    print OUT "File generated by mm_pbsa.pl\n";
    print OUT " &data\n";
    ### print OUT "  ntx    = 0,\n"; # Binary format is no longer supported
    print OUT "  ntrun  = 1,       nvect  = 0,\n";
    print OUT "  drms   = $nmdrms,\n";
    print OUT "  dielc  = 4.0,     idiel  = 0,\n";
    print OUT "  scnb   = 2.0,     scee   = 1.2\n";
    print OUT " &end\n";
    print OUT "END\n";
    close(OUT);
  }
}

sub create_sander_input(){
##########################

  # Parameters: \%GEN,\%DEC,\%SAN,\%DEL,\%GBO
  my $r_gen = shift;
  my $r_dec = shift;
  my $r_san = shift;
  my $r_del = shift;
  my $r_gbo = shift;

  print "    Sander input\n";

  my $mm = $r_gen->{"MM"} > 0 ? " MM" : "";
  my ($gb, $igb, $eedmeth, $dielc, $saltcon, $extdiel, $intdiel, $offset, $ivcap, $cutcap);
  if($r_gen->{"GB"} > 0){
    # Do GB calc
    $gb = " GB";
    $igb = $r_gbo->{"IGB"};
    if($igb == 1 || $igb == 2){
      $offset = 0.09;
    }
    else{
      $offset = 0.00;
    }
    $eedmeth = 1;
    $saltcon = $r_gbo->{"SALTCON"};
    $extdiel = $r_gbo->{"EXTDIEL"};
    $intdiel = $r_gbo->{"INTDIEL"};
    $dielc   = $r_san->{"DIELC"};
  }
  else{
    # Calc ele with const dielectric
    $gb = "";
    $igb = 0;
    $offset = 0.09;
    $eedmeth = 4;
    $saltcon = 0.0;
    $extdiel = 78.5;
    $intdiel = 1.0;
    $dielc   = $r_san->{"DIELC"};
  }
  my ($ms, $ims);
  if($r_gen->{"GB"} > 0 && $r_gen->{"MS"} == 0){
    # Calc SAS with sander
    $ms = " MS";
    $ims = $r_gbo->{"GBSA"};
  }
  else{
    # External SAS calc
    $ms = "";
    $ims = 0;
  }
  my $idecomp;
  if($r_gen->{"DC"} > 0){
    $idecomp = $r_dec->{"DCTYPE"};
  }
  else{
    $idecomp = 0;
  }

  my $suf_list = "";
  $suf_list .= "_com " if($r_gen->{"COMPLEX"});
  $suf_list .= "_rec " if($r_gen->{"RECEPTOR"});
  $suf_list .= "_lig " if($r_gen->{"LIGAND"});
  my $suf;
  foreach $suf (split/ +/,$suf_list){
    my $name = "sander" . $suf . ".in";
    open(OUT,">$name") || die("    $name not opened\n");
    print OUT "File generated by mm_pbsa.pl. Using " . $mm . $gb . $ms . "\n";
    print OUT " &cntrl\n";
    print OUT "  ntf    = 1,       ntb    = 0,       dielc  = ${dielc},\n";
    print OUT "  idecomp= ${idecomp},\n";
    print OUT "  igb    = ${igb},       saltcon= ${saltcon},\n";
    print OUT "  offset = ${offset},\n";
    print OUT "  intdiel= ${intdiel},     extdiel= ${extdiel},\n";
    # Surften is set to 1.0, so SAS is calc itself
    #   -> calc of E(SAS) is done in statistics module
    print OUT "  gbsa   = ${ims},       surften= 1.0,\n";
    print OUT "  cut    = 999.0,   nsnb   = 99999,\n";
    print OUT "  scnb   = 2.0,     scee   = 1.2,\n";
    print OUT "\n";
    print OUT "  imin   = 1,       maxcyc = 1,       ncyc   = 0,\n";
    print OUT " &end\n";
    if($eedmeth != 1){
      print OUT " &ewald\n";
      print OUT "  eedmeth  = ${eedmeth},\n";
      print OUT " &end\n";
    }

    if($r_gen->{"DC"}){
      if($suf eq "_com"){
        print OUT "Residues considered as REC\n";
        my $tmp = $r_dec->{"COMREC"};
        my @tmp = split(' ', $tmp);
        print OUT "RRES";
        my $limit = 54;
        my $length = 4;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "RRES", " ", $temp;
            $length = 5 + length $temp;
          }
        }
        print OUT "\n";
        print OUT "END\n";
        print OUT "Residues considered as LIG\n";
        $tmp = $r_dec->{"COMLIG"};
        @tmp = split(' ', $tmp);
        print OUT "LRES";
        $limit = 54;
        $length = 4;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "LRES", " ", $temp;
            $length = 5 + length $temp;
          }
        }
        print OUT "\n";
        print OUT "END\n";
        print OUT "Residues to print\n";
        $tmp = $r_dec->{"COMPRI"};
        @tmp = split(' ', $tmp);
        print OUT "RES";
        $limit = 54;
        $length = 3;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "RES", " ", $temp;
            $length = 4 + length $temp;
          }
        }
        print OUT "\n";
        print OUT "END\n";
      }
      elsif($suf eq "_rec"){
        print OUT "Residues considered as REC\n";
        my $tmp = $r_dec->{"RECRES"};
        my @tmp = split(' ', $tmp);
        print OUT "RRES";
        my $limit = 54;
        my $length = 4;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "RRES", " ", $temp;
            $length = 5 + length $temp;
          }
        }
        print OUT "\n";
        print OUT "END\n";
        print OUT "Residues to print\n";
        $tmp = $r_dec->{"RECPRI"};
        @tmp = split(' ', $tmp);
        print OUT "RES";
        $limit = 54;
        $length = 3;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "RES", " ", $temp;
            $length = 4 + length $temp;
          }
        }
        print OUT "\n";
        print OUT "END\n";
      }
      elsif($suf eq "_lig"){
        print OUT "Residues considered as LIG\n";
        my $tmp = $r_dec->{"LIGRES"};
        my @tmp = split(' ', $tmp);
        print OUT "LRES";
        my $limit = 54;
        my $length = 4;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "LRES", " ", $temp;
            $length = 5 + length $temp;
          }
        }
        print OUT "\n";
        print OUT "END\n";
        print OUT "Residues to print\n";
        $tmp = $r_dec->{"LIGRES"};
        @tmp = split(' ', $tmp);
        print OUT "RES";
        $limit = 54;
        $length = 3;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "RES", " ", $temp;
            $length = 4 + length $temp;
          }
        }
        print OUT "\nEND\n";
      }
      print OUT "END\n";
    }
    close(OUT);
  }
}

sub create_delphi_input(){
##########################

  # If REFE == 0, the reference state is assumed to be INDI/INDI.
  #   Hence, only one delphi calc with EXDI/INDI is performed,
  #     and the result taken is the reaction field energy. 
  #   Consider scaling of coulombic interactions as well.
  #   If ISTRNG > 0.0, it is necessary to calc two runs:
  #     1) ISTRNG = 0.0, getting solvation and total energy
  #     2) ISTRNG > 0.0, getting total energy
  #   Salt contribution can then be calc for LPB as:
  #     saltcontrib = total(ISTRNG > 0.0) - total(ISTRNG = 0.0).
  #   The solvation free energy for ISTRNG > 0.0 is obtained as:
  #     solvation(ISTRNG > 0.0) = solvation(ISTRNG = 0.0) + saltcontrib.
  #   Direct calc as total(ISTRNG > 0.0, EXDI > 1.0) - total(ISTRNG = 0.0, EXDI = 1.0)
  #     works for small molecules, but yields "weird" (why?) values in the latter
  #     case for big molecules.
  #   For NLPB, see work of Tsui & Case, J. Phys. Chem. B 2001, 105, 11314
  #     and salt_contrib and process_phi scripts in delphi directory.
  #
  # If REFE > 0, two delphi calc. are performed, with EXDI/INDI and
  #   1.0/INDI. The result is taken from the difference of the
  #     total energies.
  #   Hence, coulombic interactions are not taken from sander.
  #   ISTRNG > 0.0 should be straight-forward in this case.
  #

  # Parameters: \%DEL
  my $r_del = shift;
  
  print "    Delphi input\n";

  my $focruns = $r_del->{"FOCUS"} + 1;
  my $i;
  for($i = 0; $i < $focruns; $i++){

    my $multruns = 1;
    if((exists $r_del->{"ISTRNG"} && $r_del->{"ISTRNG"} > 0.0) || $r_del->{"REFE"} > 0){
      $multruns = 2;
    }    
    my $j;
    for($j = 0; $j < $multruns; $j++){

      my $name = "delphi.prm" . "." . $i . "." . $j;
      open(OUT,">$name") || die("$name not opened\n");

      my $key;
      foreach $key (keys %$r_del){
        my $prkey = lc $key;
        if($key eq "SCALE" || $key eq "PERFIL"){
          print OUT $prkey . "=" . $r_del->{$key}->[$i] . "\n";
        }
        elsif($key eq "ISTRNG"){
          if($r_del->{"REFE"} == 0){
            if($j == 0){ 
              # Zero salt in first run to get solv and total
              print OUT $prkey . "=0.0\n"; 
            }
            else{
              # ISTRNG > 0.0 
              print OUT $prkey . "=" . $r_del->{$key} . "\n";
            }
          }
          else{
            print OUT $prkey . "=" . $r_del->{$key} . "\n";
          }
        }
        elsif($key eq "EXDI" || $key eq "INDI"){
          if($r_del->{"REFE"} == 0){
            print OUT $prkey . "=" . $r_del->{$key} . "\n";
          }
          else{
            if($j == 0 && $key eq "EXDI"){
              # Calc 1.0 / INDI
              print OUT $prkey . "=" . "1.0\n";
            }
            else{
              print OUT $prkey . "=" . $r_del->{$key} . "\n";
            }
          }
        }
        elsif($key ne "REFE" &&
              $key ne "CHARGE" && 
              $key ne "SIZE" && 
              $key ne "FOCUS" &&
              $key ne "SURFTEN" &&
              $key ne "SURFOFF"){
          # Write all other parameters
          print OUT $prkey . "=" . $r_del->{$key} . "\n";
        }
      }

      # Boundary conditions
      if($i == 0){
	# Sum of Debye-Hueckel potentials
	print OUT "BNDCON=4\n";
      }
      else{
	# Focussing for runs > 0
	print OUT "BNDCON=3\n"; 
      }

      # Charge, size and pdb structure
      print OUT "in(crg,file=\"" . $r_del->{"CHARGE"} . "\")\n";
      print OUT "in(siz,file=\"" . $r_del->{"SIZE"} . "\")\n";
      print OUT "in(pdb)\n";

      # Read and write potential maps if focussing
      if($r_del->{"FOCUS"}){
        if($i > 0){
          print OUT "in(phi,unit=14)\n";
        }
        if($i < $focruns - 1){
          print OUT "out(phi,unit=14)\n";
        }
      }    

      # Output result in last run
      if($i == $focruns - 1){
        print OUT "energy(s,g)\n";
      }

      close(OUT);
    }
  }
  return;
}

sub create_pbsa_input(){
##########################

  # Parameters: \%GEN,\%DEC,\%SAN,\%DEL
  my $r_gen = shift;
  my $r_dec = shift;
  my $r_san = shift;
  my $r_del = shift;

  print "    PBSA input\n";

  my ($igb, $dielc, $idecomp, $epsin, $epsout, $istrng, $radiopt, $sprob, $space, $maxitn,
      $npopt, $cavity_surften, $cavity_offset,
      $ivcap, $cutcap, $xcap, $ycap, $zcap);

  # Do PB calculation
  $igb = 10;
  $dielc = $r_san->{"DIELC"};
  if($r_gen->{"DC"} > 0){
    $idecomp = $r_dec->{"DCTYPE"};
  }
  else{
    $idecomp = 0;
  }
  $epsin = $r_del->{"INDI"};
  $epsout = $r_del->{"EXDI"};
  if(exists $r_del->{"ISTRNG"}){
    $istrng = $r_del->{"ISTRNG"};
  }
  else{
    $istrng = 0.0;
  }
  $radiopt = $r_del->{"RADIOPT"};
  $sprob = $r_del->{"PRBRAD"};
  $space = 1.0 / $r_del->{"SCALE"}->[0];
  $maxitn = $r_del->{"LINIT"};
  $npopt = $r_del->{"NPOPT"};
  $cavity_surften = $r_del->{"SURFTEN"};
  $cavity_offset = $r_del->{"SURFOFF"};
  if(exists $r_del->{"IVCAP"} && $r_del->{"IVCAP"} == 1){
    $ivcap = 1;
    $cutcap = $r_del->{"CUTCAP"};
    $xcap = $r_del->{"XCAP"};
    $ycap = $r_del->{"YCAP"};
    $zcap = $r_del->{"ZCAP"};
  }
  elsif(exists $r_del->{"IVCAP"} && $r_del->{"IVCAP"} == 5){
    $ivcap = 5;
    $cutcap = $r_del->{"CUTCAP"};
    $xcap = 0.0;
    $ycap = 0.0;
    $zcap = 0.0;
  }
  else{
    $ivcap = 0;
    $cutcap = -1.0;
    $xcap = 0.0;
    $ycap = 0.0;
    $zcap = 0.0;
  }

  my $suf_list = "";
  $suf_list .= "_com " if($r_gen->{"COMPLEX"});
  $suf_list .= "_rec " if($r_gen->{"RECEPTOR"});
  $suf_list .= "_lig " if($r_gen->{"LIGAND"});
  my $suf;
  foreach $suf (split/ +/,$suf_list){
    my $name = "pbsa" . $suf . ".in";
    open(OUT,">$name") || die("    $name not opened\n");
    print OUT "File generated by mm_pbsa.pl. Using PB\n";
    print OUT " &cntrl\n";
    print OUT "  ntf    = 1,       ntb    = 0,\n";
    print OUT "  igb    = ${igb},      dielc  = ${dielc},\n";
    print OUT "  cut    = 999.0,   nsnb   = 99999,\n";
    print OUT "  scnb   = 2.0,     scee   = 1.2,\n";
    print OUT "  imin   = 1,       maxcyc = 0,       ntmin  = 2,\n";
    print OUT "  ivcap  = ${ivcap},       cutcap = ${cutcap},\n";
    print OUT "  xcap   = ${xcap}, ycap   = ${ycap}, zcap   = ${zcap},\n";
    print OUT "  idecomp= ${idecomp}\n";
    print OUT " &end\n";
    # Here, only REFE == 0 calculations are allowed, whereby the reference state
    #   is assumed to be INDI/INDI.
    #   Hence, only one pbsa calc with EXDI/INDI is performed,
    #     and the result taken is the reaction field energy 
    #     (given as "ELE" in the pbsa output).
    #   Consider scaling of coulombic interactions as well (by DIELC) if using INDI != 1.
    # If ISTRNG > 0.0, only one calculation is done here (compared to delphi).
    #   For NLPB, see work of Tsui & Case, J. Phys. Chem. B 2001, 105, 11314
    #     and salt_contrib and process_phi scripts in delphi directory.
    # Here, a PBTEMP of 300K is used to compute the Boltzmann factor for the salt effect.
    #   If a different temperature is specified here, also correct the factor for
    #     kcal/mol -> kT conversion in the subroutine calc_pbsa.
    print OUT " &pb\n";
    print OUT "  epsin  = ${epsin},     epsout  = ${epsout},\n";
    print OUT "  istrng = ${istrng},       radiopt = ${radiopt},\n";
    print OUT "  sprob  = ${sprob},     space   = ${space},\n";
    print OUT "  maxitn = ${maxitn},    npopt   = ${npopt},    dbfopt  = 1,\n";
    print OUT "  cavity_surften = ${cavity_surften},    cavity_offset   = ${cavity_offset},\n";
    print OUT "  npbverb= 1,\n";
    print OUT "  pbtemp = 300.0,\n";
    print OUT "  fillratio=4.0\n";
    print OUT " &end\n";

    if($r_gen->{"DC"}){
      if($suf eq "_com"){
        print OUT "Residues considered as REC\n";
        my $tmp = $r_dec->{"COMREC"};
        my @tmp = split(' ', $tmp);
        print OUT "RRES";
        my $limit = 54;
        my $length = 4;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "RRES", " ", $temp;
            $length = 5 + length $temp;
          }
        }
        print OUT "\n";
        print OUT "END\n";
        print OUT "Residues considered as LIG\n";
        $tmp = $r_dec->{"COMLIG"};
        @tmp = split(' ', $tmp);
        print OUT "LRES";
        $limit = 54;
        $length = 4;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "LRES", " ", $temp;
            $length = 5 + length $temp;
          }
        }
        print OUT "\n";
        print OUT "END\n";
        print OUT "Residues to print\n";
        $tmp = $r_dec->{"COMPRI"};
        @tmp = split(' ', $tmp);
        print OUT "RES";
        $limit = 54;
        $length = 3;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "RES", " ", $temp;
            $length = 4 + length $temp;
          }
        }
        print OUT "\n";
        print OUT "END\n";
      }
      elsif($suf eq "_rec"){
        print OUT "Residues considered as REC\n";
        my $tmp = $r_dec->{"RECRES"};
        my @tmp = split(' ', $tmp);
        print OUT "RRES";
        my $limit = 54;
        my $length = 4;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "RRES", " ", $temp;
            $length = 5 + length $temp;
          }
        }
        print OUT "\n";
        print OUT "END\n";
        print OUT "Residues to print\n";
        $tmp = $r_dec->{"RECPRI"};
        @tmp = split(' ', $tmp);
        print OUT "RES";
        $limit = 54;
        $length = 3;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "RES", " ", $temp;
            $length = 4 + length $temp;
          }
        }
        print OUT "\n";
        print OUT "END\n";
      }
      elsif($suf eq "_lig"){
        print OUT "Residues considered as LIG\n";
        my $tmp = $r_dec->{"LIGRES"};
        my @tmp = split(' ', $tmp);
        print OUT "LRES";
        my $limit = 54;
        my $length = 4;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "LRES", " ", $temp;
            $length = 5 + length $temp;
          }
        }
        print OUT "\n";
        print OUT "END\n";
        print OUT "Residues to print\n";
        $tmp = $r_dec->{"LIGRES"};
        @tmp = split(' ', $tmp);
        print OUT "RES";
        $limit = 54;
        $length = 3;
        foreach my $temp (@tmp){
          $temp =~ s/-/ /g;
          $length += 1 + length $temp;
          if($length < $limit){
            print OUT " ", $temp;
          }
          else{
            print OUT "\n", "RES", " ", $temp;
            $length = 4 + length $temp;
          }
        }
        print OUT "\n";
        print OUT "END\n";
      }
      print OUT "END\n";
    }
    close(OUT);
  }
  return;
}

sub create_iapbs_input() {
#########################
  # Parameters: \%GEN,\%SAN,\%DEL
  my $r_gen = shift;
  my $r_san = shift;
  my $r_del = shift;

  print "    iAPBS input\n";

  my ($igb, $dielc, $epsin, $epsout, $istrng, $radiopt, 
      $space, $maxitn, $npopt,
      $srad, $bcfl, $srfm, $chgm, $swin, $gamma);

  # Do iAPBS calculation
  $igb = 6;
  $dielc = $r_san->{"DIELC"};
  $epsin = $r_del->{"INDI"};
  $epsout = $r_del->{"EXDI"};
  my $pbions = "";
  if(exists $r_del->{"ISTRNG"} && $r_del->{"ISTRNG"} > 0.0){
      my $ionc = $r_del->{"ISTRNG"};
      $pbions =  "  nion = 2, \n";
      $pbions .= "  ionq  = 1.0, -1.0,\n";
      $pbions .= "  ionrr = 2.0, 2.0,\n";
      $pbions .= "  ionc  = $ionc, $ionc,\n";
  } else {
      $pbions =  "  nion = 0,\n";
  }
  $radiopt = $r_del->{"RADIOPT"};
  $srad = $r_del->{"PRBRAD"};
  $space = 1.0 / $r_del->{"SCALE"}->[0];
  $bcfl = $r_del->{"BCFL"};
  $srfm = $r_del->{"SRFM"};
  $chgm = $r_del->{"CHGM"};
  $swin = $r_del->{"SWIN"};
  $gamma = $r_del->{"GAMMA"};
 # ($dimex, $dimey, $dimez) = split(/,/,$r_del->{"DIME"});
 # ($fglenx, $fgleny, $fglenz) = split(/,/,$r_del->{"FGLEN"});
 # ($cglenx, $cgleny, $cglenz) = split(/,/,$r_del->{"CGLEN"});

  my $name = "iapbs.in";
  open(OUT,">$name") || die("    $name not opened\n");
  print OUT "File generated by mm_pbsa.pl. Using iAPBS\n";
  print OUT " &cntrl\n";
  print OUT "  ntf    = 1,       ntb    = 0,\n";
  print OUT "  igb    = ${igb},       dielc  = ${dielc},\n";
  print OUT "  cut    = 999.0,   nsnb   = 99999,\n";
  print OUT "  scnb   = 2.0,     scee   = 1.2,\n";
  print OUT "  imin   = 1,       maxcyc = 0,       ntmin  = 2,\n";
  print OUT " &end\n";
  # Here, only REFE == 0 calculations are allowed
  #
  print OUT " &apbs\n";
  print OUT "  apbs_print=0, apbs_debug=0, calc_type = 0,\n";
  print OUT "  grid = ${space}, ${space}, ${space},\n";
  if ($radiopt > 0 ){
      print OUT "  radiopt = ${radiopt}, pqr='pqr', \n";
  }
  print OUT "  sdie = ${epsout}, pdie = ${epsin}, srad=${srad}, \n";
  print OUT "  cmeth=1, bcfl=${bcfl}, srfm=${srfm}, chgm=${chgm},\n";
  print OUT $pbions;
  print OUT "  swin=${swin}, gamma=${gamma}, \n";
  print OUT "  calcenergy=1, calcforce=0, calcnpenergy=1 \n";
  print OUT " &end\n";

  close(OUT);
  return;
}

sub create_makecrd_input(){
###########################

  # Parameters: \%GEN,\%MAK,\$TRA
  my $r_gen = shift;
  my $r_mak = shift;
  my $r_tra = shift;

  print "    make_crd input\n";

  open(OUT,">make_crd.in") || die("make_crd.in not opened\n");
  
  if($r_gen->{"AS"}){
    my $groups = $r_mak->{"NUMBER_MUTANT_GROUPS"};
    print OUT "NUMBER_MUTANT_GROUPS=" . $groups . "\n";

    my $i;
    for($i = 0; $i < $groups; $i++){
      if(defined $r_mak->{"MUTANT_ATOM1"}->[$i]){
        print OUT "MUTANT_ATOM1=" . $r_mak->{"MUTANT_ATOM1"}->[$i] . "\n";
      }
      else{
        die("Something wrong with MUTANT_ATOM1 definitions\n");
      }

      if(defined $r_mak->{"MUTANT_ATOM2"}->[$i]){
        print OUT "MUTANT_ATOM2=" . $r_mak->{"MUTANT_ATOM2"}->[$i] . "\n";
      }
      else{
        die("Something wrong with MUTANT_ATOM2 definitions\n");
      }

      if(defined $r_mak->{"MUTANT_KEEP"}->[$i]){
        print OUT "MUTANT_KEEP=" . $r_mak->{"MUTANT_KEEP"}->[$i] . "\n";
      }
      else{
        die("Something wrong with MUTANT_KEEP definitions\n");
      }

      if(defined $r_mak->{"MUTANT_REFERENCE"}->[$i]){
        print OUT "MUTANT_REFERENCE=" . $r_mak->{"MUTANT_REFERENCE"}->[$i] . "\n";
      }
      else{
        die("Something wrong with MUTANT_REFERENCE definitions\n");
      }
    } # end for
  } # end if

  my $groups = $r_mak->{"NUMBER_LIG_GROUPS"};
  print OUT "NUMBER_LIG_GROUPS=" . $groups . "\n";
  my $i;
  for($i = 0; $i < $groups; $i++){
    if(defined $r_mak->{"LSTART"}->[$i]){
      print OUT "LSTART=" . $r_mak->{"LSTART"}->[$i] . "\n";
    }
    else{
     die("Something wrong with LSTART definitions\n");
    }
     if(defined $r_mak->{"LSTOP"}->[$i]){
      print OUT "LSTOP=" . $r_mak->{"LSTOP"}->[$i] . "\n";
    }
    else{
      die("Something wrong with LSTOP definitions\n");
    }
  } # end for
        
  $groups = $r_mak->{"NUMBER_REC_GROUPS"};
  print OUT "NUMBER_REC_GROUPS=" . $groups . "\n";
  for($i = 0; $i < $groups; $i++){
    if(defined $r_mak->{"RSTART"}->[$i]){
      print OUT "RSTART=" . $r_mak->{"RSTART"}->[$i] . "\n";
    }
    else{
      die("Something wrong with RSTART definitions\n");
    }
    if(defined $r_mak->{"RSTOP"}->[$i]){
      print OUT "RSTOP=" . $r_mak->{"RSTOP"}->[$i] . "\n";
    }
    else{
      die("Something wrong with RSTOP definitions\n");
    }
  } # end for

  my $key;
  foreach $key (@GC_PAR){
    if($key ne "NUMBER_LIG_GROUPS" && $key ne "NUMBER_REC_GROUPS"){
      print OUT $key . "=" . $r_mak->{$key} . "\n";
    }
  }

  print OUT "PREFIX=" . $r_gen->{"PATH"} . $r_gen->{"PREFIX"} . "\n";
  if($r_gen->{"COMPLEX"}){
    print OUT "SUFFIX_C=_com\n";
  }
  else{
    print OUT "SUFFIX_C=\n";
  }
  if($r_gen->{"RECEPTOR"}){
    print OUT "SUFFIX_R=_rec\n";
  }
  else{
    print OUT "SUFFIX_R=\n";
  }
  if($r_gen->{"LIGAND"}){
    print OUT "SUFFIX_L=_lig\n";
  }
  else{
    print OUT "SUFFIX_L=\n";
  }

  close(OUT);
  return;
}

1; # Necessary for package function