#!/usr/local/bin/perl -w

# Author: Jan Labanowski, done on 97.03.23
#
# this utility converts orbital basis sets from DGauss (compact) format
# deMon format which is easier for humans to read. It has to be used as
# a filter, e.g.: 
#     cat DGauss.bas | DGauss2deMon-o.pl > deMon.bas
# Remember that first line of this file has to point to the location of perl
# interpreter, and this very file has to have execute permission  bit(s) set.
# The DGauss.bas is a file containing a single definition of a basis
# set for some element, and deMon.bas is a file which will contain
# the same basis in the deMon format.
# If you find bugs, please bug Jan Labanowski, jkl@ccl.net

%elements = (
 'H',  'HYDROGEN',
 'HE', 'HELIUM',
 'LI', 'LITHIUM',
 'BE', 'BERYLLIUM',
 'B',  'BORON',
 'C',  'CARBON',
 'N',  'NITROGEN',
 'O',  'OXYGEN',
 'F',  'FLURINE',
 'NE', 'NEON',
 'NA', 'SODIUM',
 'MG', 'MAGNESIUM',
 'AL', 'ALUMINIUM',
 'SI', 'SILICON',
 'P',  'PHOSPHORUS',
 'S',  'SULFUR',
 'CL', 'CHLORINE',
 'AR', 'ARGON',
 'K',  'POTASSIUM',
 'CA', 'CALCIUM',
 'SC', 'SCANDIUM',
 'TI', 'TITANIUM',
 'V',  'VANADIUM',
 'CR', 'CHROMIUM',
 'MN', 'MANGANESE',
 'FE', 'IRON',
 'CO', 'COLBALT',
 'NI', 'NICKEL',
 'CU', 'COPPER',
 'ZN', 'ZINC',
 'GA', 'GALLIUM',
 'GE', 'GERMANIUM',
 'AS', 'ARSENIC',
 'SE', 'SELENIUM',
 'BR', 'BROMIUM',
 'KR', 'KRYPTON',
 'RB', 'RUBIDIUM',
 'SR', 'STRONTIUM',
 'Y',  'YTtRIUM',
 'ZR', 'ZIRCONIUM',
 'NB', 'NOBIUM',
 'MO', 'MOLYBDENUM',
 'TC', 'TECHNICIUM',
 'RU', 'RUTHENIUM',
 'RH', 'RHODIUM',
 'PD', 'PALLADIUM',
 'AG', 'SILVER',
 'CD', 'CADIUM',
 'IN', 'INDIUM',
 'SN', 'TIN',
 'SB', 'ANTYMONIUM',
 'TE', 'TELLIUM',
 'I',  'IODINE',
 'XE', 'XENON');

$line = &getline();
if($line !~ /^(\S+)\s+(\S+)/) {
  die "This line should have a format\nEL  BASISTYPE\n";
  }
$symbol = $1;
$basis = $2;
$line = &getline();
$line =~ s/^\s+//;
@fields = split(/\s+/, $line);
$n_types = 0;
for($i = 0; $i <= $#fields; $i++) {
  if($fields[$i] =~ /(\d+)/) {
    $n_types++;
    $types[$n_types] = int($1);
    }
  }
if($n_types <= 0) {
  die "The line with number of contractions for s,p,d not given:\n$line\n";
  }

$n_contr = 0;
$n_exp = 0;
$n_coefs = 0;
for ($n = 1; $n <= $n_types; $n++) {
  for ($j = 1; $j <= $types[$n]; $j++) {
    $n_contr++;
    $contraction_starts[$n_contr] = $n_exp + 1;
    $cont = 1;
    while ($cont) {
      $line = &getline();
      if($line =~ /^\s*\S+=([\d]+\.[\d]+)\s*$/) {
        $scale = $1 + 0.0;
        $line = &getline();
        }
      else {
        $scale = 1.0;
        }
      @fields = split(/\s+/, $line);
      for($i = 0; $i <= $#fields; $i++) {
        if($fields[$i] =~ /:$/) {  # skip comments
          next;
          }
        if($fields[$i] =~ /^[-+E.0-9]+$/) {
          $n_exp++;
          $exponents[$n_exp] = $fields[$i]*$scale*$scale;
          }
        else {
          die "Non-number on the exponents line:\n$line\n";
          }
        }
      if(($#fields == 0) && (($n_coefs+1) == $contraction_starts[$n_contr])) {
#      if($#fields == 0) {
        $n_coefs++;
        $coefs[$n_coefs] = 1.0;
        $cont = 0;
        }
      else {
        $line = &getline();    # coeficients line
        @fields = split(/\s+/, $line);
        $cont = 0;
        for($i = 0; $i <= $#fields; $i++) {
          if($fields[$i] =~ /:$/) {  # skip comments
            next;
            }
          if($fields[$i] =~ /^[-+E.0-9]+$/) {
            $n_coefs++;
            $coefs[$n_coefs] = $fields[$i] + 0.0;
            }
          elsif ($i == $#fields) {
            $cont = 1;
            }
          else {
            die "Non-number on the coefficients line:\n$line\n";
            }
          }
        }
      }
    
    if($n_coefs != $n_exp) {
      die "No. of exponents not equal to no. of coefficients:\n$line\n";
      }
    $contraction_ends[$n_contr] = $n_exp;
    }
  }

$n_contr = 0;
$label = 'O-' . $elements{$symbol} . ' (';
for ($n = 1; $n <= $n_types; $n++) {
  for ($j = 1; $j <= $types[$n]; $j++) {
    $n_contr++;
    $n_gauss =$contraction_ends[$n_contr] - $contraction_starts[$n_contr] + 1;
    $label .= sprintf("%1d", $n_gauss);
    }
  if($n < $n_types) {
    $label .= '/';
    }
  }
$label .= ')';
print STDOUT "# $basis\n";
print STDOUT "$label\n";
for ($n = $n_types+1; $n <= 3; $n++) {   # if no p's and/or d's
  $types[$n] = 0;
  }

for ($n = 1; $n <= 3; $n++) {  # print number of s, p, & d contractions
  printf STDOUT "%5d", $types[$n];
  }
print STDOUT "\n";

$n_contr = 0;
for ($n = 1; $n <= $n_types; $n++) {
  for ($j = 1; $j <= $types[$n]; $j++) {
    $n_contr++;
    $n_gauss = $contraction_ends[$n_contr] - $contraction_starts[$n_contr] + 1;
    printf STDOUT "%5d\n", $n_gauss;        
    for ($i = $contraction_starts[$n_contr];
         $i <= $contraction_ends[$n_contr]; $i++) { # print exps and coefs
      printf STDOUT "%20.10f%20.10f\n", $exponents[$i], $coefs[$i];
      }
    }
  }



# =====================================================
sub getline {
  local ($inpl);
  local (@tokens);
  local ($nt);
   while($inpl = <STDIN>) {
     if(($inpl !~ /\S/) || ($inpl =~ /^\s*:/)) {
       next;
       }
     else {
       chop($inpl);
       $inpl =~ s/^\s+//;
       $inpl =~ s/\s+$//;
       @tokens = split(/\s+/, $inpl);
       $inpl = "";
       for($nt = 0; $nt <= $#tokens; $nt++) {
         if($tokens[$nt] =~ /:$/) {
           next;
           }
         else {
           if($inpl ne "") {
             $inpl = $inpl . " ";
             }
           $inpl = $inpl . $tokens[$nt];
           }
         }
       if($inpl ne "") {
         last;
         }
       else {
         next;
         }
       }
     }
   return($inpl);
   }