#!/usr/bin/perl -w
use strict;

my $description = <<'END_DESC';
# Reads an ISIM input file (e.g. "run.in"), along with the _EXCESS files
# associated with the input file.
# Writes an xml file of ion conditions to standard output.
# Converts soft atoms to hard sphere, to keep ISIM from being so slow.
END_DESC

my $program_name = $0;
my $usage = "Usage: $program_name <isim_parameter_input_file1> [<isim_parameter_input_file2] ...\n\n$description\n";

die $usage unless $#ARGV >= 0;

my $float_regexp = '[-+0-9\.]+';

print <<'END_HEAD';
<?xml version="1.0"?>
<!DOCTYPE salt_condition_range_list SYSTEM "salt_condition_range_list.dtd">
<salt_condition_range_list>
END_HEAD

# Open each file mentioned on the command line
foreach my $isim_parameter_input (@ARGV) {
  open FH, "<$isim_parameter_input" or die "Unable to open file $isim_parameter_input";

  print STDERR "Parsing file $isim_parameter_input...\n";

  # 1) Parse temperature, solvent permittivity, and ion names
  my ($temperature, $solvent_permittivity);
  my $ion_surface_model = "HARD_SPHERE";
  my @ions = ();

  my $is_in_ion_section = 0;
  while (<FH>) {
    if ($_ !~ /\S/) {next;} # blank lines

    elsif ($_ =~ /^\s*#ions\s*$/i) {
      $is_in_ion_section = 1;
    }

    elsif ($_ =~ /^\s*#endions\s*$/i) {
      $is_in_ion_section = 0;
    }

    # temperature
    elsif ($_ =~ /^\s*temperature\s+($float_regexp)\s*$/i) {
      $temperature = $1;
    }

    # solvent permittivity
    elsif ($_ =~ /^\s*medium_permittivity\s+($float_regexp)\s*$/i) {
      $solvent_permittivity = $1;
    }

    # ion surface model
    elsif ($_ =~ /^\s*model\s+(L|H)\s*$/i) {
      if ($1 eq "H") {$ion_surface_model = "HARD_SPHERE";}
      else {$ion_surface_model = "LENNARD_JONES";}
      $is_in_ion_section = 1;
    }

    # number of ion types
    elsif ($_ =~ /^\s*types\s+(\d+)\s*$/i) {
      $is_in_ion_section = 1;
    }

    elsif ($is_in_ion_section) {
      # if we get here, we expect an ion parameter line

      # Format depends upon whether its hard sphere or soft
      my ($ion_name, $concentration, $charge, $radius, $well_depth);
      if ($ion_surface_model eq "HARD_SPHERE") {
        die $_ unless $_ =~ /^
          \s*(\S+)  # ion name
          \s+($float_regexp) # concentration
          \s+($float_regexp) # charge
          \s+($float_regexp) # radius
          \s*$/x;
        $ion_name = $1;
        $concentration = $2;
        $charge = $3;
        $radius = $4;
      }
      else { # Lennard Jones
        die unless $_ =~ /^
          \s*(\S+)  # ion name
          \s+($float_regexp) # concentration
          \s+($float_regexp) # charge
          \s+($float_regexp) # radius
          \s+($float_regexp) # well depth
          \s*$/x;
        $ion_name = $1;
        $concentration = $2;
        $charge = $3;
        $radius = $4;
        $well_depth = $5;
      }

      # Populate ion data
      my %ion = ();
      $ion{NAME} = $ion_name;
      $ion{CONCENTRATION} = $concentration;
      $ion{CHARGE} = $charge;
      $ion{RADIUS} = $radius;
      $ion{WELL_DEPTH} = $well_depth if defined $well_depth;

      push @ions, \%ion;
    }
  }

  die "Unable to find temperature in file $isim_parameter_input" unless defined $temperature;
  die "Unable to find solvent permittivity in file $isim_parameter_input" unless defined $solvent_permittivity;
  
  # TODO - read excess chemical potential from EXCESS files

  # Data structure for storing potential data
  my %condition_ion_excess = ();
  my %condition_ion_concentration = ();

  foreach my $ion (@ions) {
    my $excess_file_name = $ion->{NAME}."_EXCESS";

    # get directory name of input file
    my $dir_name = "";
    # Is there a directory in the file name?
    if ($isim_parameter_input =~ m!^(.*/)[^/]+$! ) {
      $dir_name = $1;
    }

    my $full_excess = $dir_name.$excess_file_name;
    open EXCESS, "<$full_excess" or die "No excess file $full_excess";

    # Read header line
    my $header_line = <EXCESS>;
    chomp $header_line;
    my $ion_count = $#ions + 1;
    my $header_regexp = '^#\s*((?:\S+\s+){'.$ion_count.'})kcal/mol\s*$';
    die $header_line unless $header_line =~ m!$header_regexp!;
    my @header_ions = split /\s+/, $1;
    # print join(", ", @header_ions) . "\n";

    # Read data lines
    while (<EXCESS>) {
      chomp $_;
      next unless $_ =~ /\S/; # ignore blank lines
      # TODO
      my $data_regexp = '^\s*((?:'.$float_regexp.'\s+){'.$ion_count.'})('.$float_regexp.')\s*$';
      die $_ unless $_ =~ m!$data_regexp!;

      # TODO

      my @concentrations = split /\s+/, $1;
      # TODO create unique key name for each set of conditions
      my %condition_ions = ();
      my $index = 0;
      foreach my $condition_ion(@header_ions) {
        $condition_ions{$header_ions[$index]} = $concentrations[$index];
        $index ++;
      }
      my $condition_key = "";
      foreach my $condition_ion (sort keys %condition_ions) {
        # concentration
        $condition_key .= sprintf("%.2f", $condition_ions{$condition_ion});

        # ion name
        $condition_key .= $condition_ion;

        # $condition_key .= " "; # separate ions with spaces
      }

      my $excess_potential = $2;

      # print "$condition_key\t$excess_potential\n";

      my $concentration = $condition_ions{$ion->{NAME}};

      $condition_ion_excess{$condition_key}{$ion->{NAME}} = $excess_potential;
      $condition_ion_concentration{$condition_key}{$ion->{NAME}} = $concentration;
    }

    close EXCESS;
  }

  # write xml fragment for this condition
  #  <salt_condition_range
  #      surface_model = "HARD_SPHERE"
  #  >
  #    <ion_species 
  #      ion_id = "Mg2plus"
  #      charge = "2.0"
  #      radius = "3.0"
  #    />
  #    <ion_species
  #      ion_id = "Cl1minus"
  #      charge = "-1.0"
  #      radius = "4.0"
  #    />
  #    <salt_condition condition_name="5mM MgCl2">

  # open salt_condition_range
  print "  <salt_condition_range\n    surface_model = \"HARD_SPHERE\"\n  >\n";
  # print "  <salt_condition_range\n    surface_model = \"$ion_surface_model\"\n  >\n";

  # One ion_species stanza for each ion type
  foreach my $ion (@ions) {
    print "    <ion_species\n";
    print "      ion_id = \"".$ion->{NAME}."\"\n";
    print "      charge = \"".$ion->{CHARGE}."\"\n";
    print "      radius = \"".$ion->{RADIUS}."\"\n";
    # print "      well_depth = \"".$ion->{WELL_DEPTH}."\"\n" if exists $ion->{WELL_DEPTH};
    print "    />\n";
  }

  # TODO - one salt_condition for each discrete excess chemical potential
  foreach my $ion_condition (sort keys %condition_ion_excess) {

    my $indent = "    ";

    # open tag
    print $indent."<salt_condition condition_name=\"$ion_condition\">\n";

    foreach my $ion (sort keys %{$condition_ion_excess{$ion_condition}}) {
      my $indent2 = $indent . "  ";

      my $excess_potential = $condition_ion_excess{$ion_condition}{$ion};
      my $concentration = $condition_ion_concentration{$ion_condition}{$ion};

      print $indent2."<ion_excess_potential ion_id=\"$ion\" concentration=\"$concentration\" excess_chemical_potential=\"$excess_potential\"/>\n";
    }

    # close tag
    print $indent."</salt_condition>\n";
  }

  # close salt_condition_range
  print "  </salt_condition_range>\n";

  close FH;
}

print "</salt_condition_range_list>\n";