#!/usr/bin/perl

# Converts GAMESS and Priroda output files
# into IGLO, CLOPPA, Priroda and Dalton input files

use Getopt::Std;

$BORHTOA = 0.52917725;

# Predefined tables...

# atom names (periodic table)
@atname = qw (XX
        H                                                  He
        Li Be                               B  C  N  O  F  Ne
        Na Mg                               Al Si P  S  Cl Ar
        K  Ca Sc Ti V  Cr Mn Fe Co Ni Cu Zn Ga Ge As Se Br Kr
        Rb Sr Y  Zr Nb Mo Tc Ru Rh Pd Ag Cd In Sn Sb Te I  Xe
        Cs Ba La
                Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu
                 Hf Ta W  Re Os Ir Pt Au Hg Tl Pb Bi Po At Rn
        Fr Ra Ac
                Th Pa U  Np Pu Am Cm Bk Cf Es Fm Md No Lr
       );

# orbital names
@orbname = qw (s p d f);

# atom basises (IGLO)
%basisdz = (
      'H',  'h03;c,1.2',
      'Li', 'h07;c,1.4',
      'Be', 'h07;c,1.4',
      'B',  'h0703;c,1.4 h0703;c,1.2',
      'C',  'h0703;c,1.4 h0703;c,1.2',
      'N',  'h0703;c,1.4 h0703;c,1.2',
      'O',  'h0703;c,1.4 h0703;c,1.2',
      'F',  'h0703;c,1.4 h0703;c,1.2',
      'Na', 'h1006;c,1.5 h1006;c,1.3',
      'Mg', 'h1006;c,1.5 h1006;c,1.3',
      'Al', 'h1006;c,1.5 h1006;c,1.3',
      'Si', 'h1006;c,1.5 h1006;c,1.3',
      'P',  'h1006;c,1.5 h1006;c,1.3',
      'S',  'h1006;c,1.5 h1006;c,1.3',
      'Cl', 'h1006;c,1.5 h1006;c,1.3'
     );
%basis1 = (
     'H',  'h05;c,1.3 0.65',
     'Li', 'h09;c,1.5 0.19,0.75,3.0',
     'Be', 'h09;c,1.5 0.058,0.18,0.88,2.2 0.5',
     'B',  'h0905;c,1.5 h0905;c,1.3 0.7',
     'C',  'h0905;c,1.5 h0905;c,1.3 1.0',
     'N',  'h0905;c,1.5 h0905;c,1.3 1.0',
     'O',  'h0905;c,1.5 h0905;c,1.3 1.0',
     'F',  'h0905;c,1.5 h0905;c,1.3 1.0',
     'Na', 'h1107;c,1.5 h1107;c,1.3 0.2,0.8',
     'Mg', 'h1107;c,1.5 h1107;c,1.3 0.25,1.0',
     'Al', 'h1107;c,1.5 h1107;c,1.3 0.3,1.2',
     'Si', 'h1107;c,1.5 h1107;c,1.3 0.35,1.4',
     'P',  'h1107;c,1.5 h1107;c,1.3 0.35,1.4',
     'S',  'h1107;c,1.5 h1107;c,1.3 0.4,1.6',
     'Cl', 'h1107;c,1.5 h1107;c,1.3 0.4,1.6'
    );
%basis2 = (
     'H',  'h05;c,1.3 0.65',
     'Li', 'h09;c,1.5 0.19,0.75,3.0',
     'Be', 'h09;c,1.5 0.058,0.18,0.88,2.2 0.5',
     'B',  'h0905;c,1.5 h0905;c,1.2 0.7',
     'C',  'h0905;c,1.5 h0905;c,1.2 1.0',
     'N',  'h0905;c,1.5 h0905;c,1.2 1.0',
     'O',  'h0905;c,1.5 h0905;c,1.2 1.0',
     'F',  'h0905;c,1.5 h0905;c,1.2 1.0',
     'Na', 'h1107;c,1.5 h1107;c,1.2 0.2,0.8',
     'Mg', 'h1107;c,1.5 h1107;c,1.2 0.25,1.0',
     'Al', 'h1107;c,1.5 h1107;c,1.2 0.3,1.2',
     'Si', 'h1107;c,1.5 h1107;c,1.2 0.35,1.4',
     'P',  'h1107;c,1.5 h1107;c,1.2 0.35,1.4',
     'S',  'h1107;c,1.5 h1107;c,1.2 0.4,1.6',
     'Cl', 'h1107;c,1.5 h1107;c,1.2 0.4,1.6'
    );
%basis3 = (
     'H', 'h06;c,1.3 0.33,1.3',
     'C', 'h1107;c,1.5 h1107;c,1.2 0.35,1.4',
     'N', 'h1107;c,1.5 h1107;c,1.2 0.35,1.4',
     'O', 'h1107;c,1.5 h1107;c,1.2 0.35,1.4',
     'F', 'h1107;c,1.5 h1107;c,1.2 0.35,1.4',
     'P', 'h1208;c,1.5 h1208;c,1.2 0.25,0.8,2.5',
     'S', 'h1208;c,1.5 h1208;c,1.2 0.25,0.8,2.5'
    );
%basis4 = (
     'H', 'h06;c,1.2 0.1,0.4,1.6 0.65',
     'C', 'h1107;c,1.4 h1107 0.2,0.8,3.2 1.0',
     'N', 'h1107;c,1.4 h1107 0.2,0.8,3.2 1.0',
     'O', 'h1107;c,1.4 h1107 0.2,0.8,3.2 1.0',
     'F', 'h1107;c,1.4 h1107 0.2,0.8,3.2 1.0',
     'P', 'h1208;c,1.4 h1208 0.15,0.6,2.4,9.6 0.5,1.5',
     'S', 'h1208;c,1.4 h1208 0.15,0.6,2.4,9.6 0.5,1.5'
    );

# fallback atom ECP and basis (IGLO) - if not found above
%basisecp = (
       'Br', 'ecp28mwb'
      );
%basisfb = (
      'Br', 'ecp28mwb;c ecp28mwb;c'
     );

# Subroutines to interpret input files...

# read GAMESS format
sub readgamess
{
  my ($mem, $mult, $pointgroup, $axisorder, $gpattern, $iatom, $atom, $charge,
      $x, $y, $z, $ilabel, $nlabels);

  # search for memory
  while (<STDIN>)
  {
    chomp;
    last if ($memory) = /^ *(\d+) +WORDS OF MEMORY AVAILABLE$/;
  }
  $memory /= 1024*1024;
  ($mem, $mult) = $opt_m =~ /^(\d*\.?\d*)([km]?)$/i;
  if ("\L$mult" eq 'm')
  {
    $mult = 1;
  }
  elsif ("\L$mult" eq 'k')
  {
    $mult = 1024;
  }
  else
  {
    $mult = 1024*1024;
  }
  $memory = $mem/$mult if $mem > 0;
  $memory = int ($memory+0.5);
  $memory = 1 if $memory < 1;
  $memory = sprintf "memory,%.0f,m", $memory;
  
  # search for title
  while (<STDIN>)
  {
    chomp;
    last if /^ *RUN TITLE$/;
  }
  <STDIN>;
  $_ = <STDIN>;
  chomp;
  ($title) = /^ *(.*?) *$/;
  
  # search for symmetry
  while (<STDIN>)
  {
    chomp;
    last if ($pointgroup) = /^ *THE POINT GROUP OF THE MOLECULE IS +(.+?) *$/;
  }
  $_ = <STDIN>;
  chomp;
  ($axisorder) = /^ *THE ORDER OF THE PRINCIPAL AXIS IS +(\d+) *$/;
  $symm = '';
  if ($pointgroup eq 'CS')
  {
    $symm = 'z';
  }
  elsif ($pointgroup eq 'CI')
  {
    $symm = 'xyz';
  }
  elsif ($pointgroup eq 'CN')
  {
    $symm = 'xy' if $axisorder == 2;
  }
  elsif ($pointgroup eq 'CNH')
  {
    $symm = 'xy,z' if $axisorder == 2;
  }
  elsif ($pointgroup eq 'CNV')
  {
    $symm = 'x,y' if $axisorder == 2;
  }
  elsif ($pointgroup eq 'DN')
  {
    $symm = 'xy,xz,yz' if $axisorder == 2;
  }
  elsif ($pointgroup eq 'DNH')
  {
    $symm = 'x,y,z' if $axisorder == 2;
  }
  $gsymm = $pointgroup;
  $gsymm =~ s/N/$axisorder/gi;
  $gsymm =~ tr/HISV/hisv/;
  $gsymm = 'C1' unless $symm;
  
  # search for charge
  while (<STDIN>)
  {
    chomp;
    last if ($icharg) = /^ *CHARGE OF MOLECULE += *(\d+) *$/;
  }
  
  # search for geometry
  $gpattern = 'COORDINATES OF ALL ATOMS ARE';
#  $gpattern = 'COORDINATES OF SYMMETRY UNIQUE ATOMS' if $symm && $opt_o ne 'c';
  $gpattern = 'COORDINATES OF SYMMETRY UNIQUE ATOMS' if $symm && $opt_o !~ /[cpm]/ ;
  $natoms = 0;
  RPATH:while (1)      # the very last coordinates will be accepted
  {
    if ($rpath >= 0)      # stationary point might be expected
    {
      $old_rpath = $rpath;
      while (<STDIN>)
      {
        chomp;
        if (/^1\s+POINT\s*\d+ ON THE REACTION PATH$/)
        {
          $IRC = 1;
          $rpath ++;      # another IRC point found
          last;
        }
        if (/^1     \*\*\*\*\* [\w ]+ LOCATED \*\*\*\*\*$/)
        {
          $rpath ++;      # another stationary point found
          last;
        }
      }
    }
    exit 0 if $old_rpath && $old_rpath == $rpath;  # no more points found
    if ($rpath == 0)  # no stationary points found - no reaction path
    {
      $rpath = -1;  # no reaction path - rewind and scan for last coords
      seek (STDIN, 0, 0);
    }

    # Для IRC сначала надо искать энергию, а затем геометрию
    if ($opt_o eq 'm' && $IRC && $rpath>0) {
      while (<STDIN>) {
        if (/TOTAL ENERGY\s*=/) {
          ($energy) = /TOTAL ENERGY\s*=\s*(-?\d*\.\d*)/;
          last;
        }
        last if (/^ *$gpattern/);
      }
    }
    # Ищем и зачитываем геометрию
    while (<STDIN>)
    {
      chomp;
      last if /^ *$gpattern/;
    }
    last unless <STDIN>;
    <STDIN>;
    $iatom = 0;
    while (<STDIN>)
    {
      chomp;
      ($atom, $charge, $x, $y, $z) = split;
      $charge = int ($charge);
      last unless $charge;
      $atom[$iatom++] = join (' ', $charge, $x, $y, $z, $atom);
      $label{$charge} = 'a';
      $natoms{$charge}++ unless $natoms;
    }
    # Для обычной оптимизации после геометрии искать энергию
    if ($opt_o eq 'm' && ! $IRC) {
      while (<STDIN>) {
        if (/ NSERCH=/) {
          ($energy) = /NSERCH=[\s\d]+ENERGY\s*=\s*(-?\d*\.\d*)/;
          last;
        }
        if (/TOTAL ENERGY\s*=/) {
          ($energy) = /TOTAL ENERGY\s*=\s*(-?\d*\.\d*)/;
          last;
        }
        last if (/^ *$gpattern/);
      }
    }

    last RPATH unless $iatom;      # no more coordinates
    $natoms = $iatom unless $natoms;
    if ($opt_o eq 'i')
    {
      @atom = sort bycharge @atom;
      $ilabel = 0;
      foreach $charge (sort bycharge keys %label)
      {
        $label{$charge} .= ++$ilabel;
      }
      $nlabels = $ilabel;
    }
    if ($rpath > 0)        # evtl create output file
    {
      if ($ARGV[1])
      {
        $point = sprintf "%03d", $rpath;
        open (STDOUT, ">$ARGV[1]$point") ||
          die "Cannot open output file $ARGV[1]$point";
      }
      geniglo    () if $opt_o eq 'i';    # generate IGLO output file
      gencloppa  () if $opt_o eq 'c';    # generate CLOPPA output file
      gendalton  () if $opt_o eq 'd';    # generate DALTON output file
      genpriroda () if $opt_o eq 'p';    # generate PRIRODA output file
      genmolden  () if $opt_o eq 'm';   # generate MOLDEN xyz file
    }
  }
}

# read Priroda format
sub readpriroda
{
  my ($mem, $mult, $iatom, $charge, $x, $y, $z, $ilabel, $nlabels, $old_rpath);

  # search for memory
  while (<STDIN>)
  {
    chomp;
    last if ($memory) = /^ *Memory: +(\d+) MB$/;
  }
  $memory /= 8;
  ($mem, $mult) = $opt_m =~ /^(\d*\.?\d*)([km]?)$/i;
  if ("\L$mult" eq 'm')
  {
    $mult = 1;
  }
  elsif ("\L$mult" eq 'k')
  {
    $mult = 1024;
  }
  else
  {
    $mult = 1024*1024;
  }
  $memory = $mem/$mult if $mem > 0;
  $memory = int ($memory+0.5);
  $memory = 1 if $memory < 1;
  $memory = sprintf "memory,%.0f,m", $memory;
  
  # search for title
  while (<STDIN>)
  {
    chomp;
    last if ($title) = /^ *formula: +(.+)$/;  # no real title in priroda out
  }
  
  # Priroda supports no symmetry
  $symm  = '';
  $gsymm = 'C1';
  
  # search for charge
  while (<STDIN>)
  {
    chomp;
    last if ($icharg) = /^ *\| +Charge of molecule +(\d+)$/;
  }
  
  # search for geometry and evtl create output file for each stationary point
  $natoms = 0;
  RPATH:while (1)  # each stationary point coordinates will be accepted
  {
    if ($rpath >= 0)      # stationary point might be expected
    {
      $old_rpath = $rpath;
      while (<STDIN>)
      {
        chomp;
        if (/^ *OPTIMIZATION CONVERGED in /)
        {
          $rpath ++;      # another stationary point found
          last;
        }
      }
    }
    exit 0 if $old_rpath && $old_rpath == $rpath;  # no more points found
    if ($rpath == 0)  # no stationary points found - no reaction path
    {
      $rpath = -1;  # no reaction path - rewind and scan for last coords
      seek (STDIN, 0, 0);
    }

    # сначала надо искать энергию, а затем геометрию
    if ($opt_o eq 'm') {
      while (<STDIN>) {
        if (/Energy ?[:=]\s*(-?\d*\.\d*)/) {
          $energy = $1;
          last;
        }
        last if (/^ *Atomic Coordinates:$/);
      }
    }
    
    while (<STDIN>)
    {
      chomp;
      last if /^ *Atomic Coordinates:$/;
    }
    $iatom = 0;
    while (<STDIN>)
    {
      chomp;
      ($charge, $x, $y, $z) = split;
      last if $charge eq '#';
      $atom[$iatom++] = join (' ', $charge, $x, $y, $z);
      $label{$charge} = 'a';
      $natoms{$charge}++ unless $natoms;
    }
    last RPATH unless $iatom;      # no more coordinates
    $natoms = $iatom unless $natoms;
    if ($opt_o eq 'i')
    {
      @atom = sort bycharge @atom;
      $ilabel = 0;
      foreach $charge (sort bycharge keys %label)
      {
        $label{$charge} .= ++$ilabel;
      }
      $nlabels = $ilabel;
    }
    if ($rpath > 0)        # evtl create output file
    {
      if ($ARGV[1])
      {
        $point = sprintf "%03d", $rpath;
        open (STDOUT, ">$ARGV[1]$point") ||
          die "Cannot open output file $ARGV[1]$point";
      }
      geniglo    () if $opt_o eq 'i';    # generate IGLO output file
      gencloppa  () if $opt_o eq 'c';    # generate CLOPPA output file
      gendalton  () if $opt_o eq 'd';    # generate DALTON output file
      genpriroda () if $opt_o eq 'p';    # generate PRIRODA output file
      genmolden  () if $opt_o eq 'm';    # generate MOLDEN xyz file
    }
  }
}

# read GAMESS IRC format from scr
sub readscrirc
{
  my ($iatom, $charge, $x, $y, $z, $ilabel, $nlabels, $old_rpath);

  # search for geometry and evtl create output file for each stationary point
  $natoms = 0;
  RPATH:while (1)  # each stationary point coordinates will be accepted
  {
    if ($rpath >= 0)      # stationary point might be expected
    {
      $old_rpath = $rpath;
      while (<STDIN>)
      {
        chomp;
        if (/^\*\*\*\*\* BEGIN IRC INFORMATION PACKET/)
        {
          $rpath ++;      # another stationary point found
          last;
        }
      }
    }
    exit 0 if $old_rpath && $old_rpath == $rpath;  # no more points found
    if ($rpath == 0)  # no stationary points found - no reaction path
    {
      $rpath = -1;  # no reaction path - rewind and scan for last coords
      seek (STDIN, 0, 0);
    }

    # сначала надо искать энергию, а затем геометрию
    if ($opt_o eq 'm') {
      while (<STDIN>) {
        if (/^POINT= *\d+ +(STOTAL= *\S+) +E= *(-?\d+\.?\d+)/) {
          $title = $1;
          $energy = $2;
          last;
        }
        last if (/^CARTESIAN COORDINATES \(BOHR\)/);
      }
    }
    
    while (<STDIN>)
    {
      chomp;
      last if /^CARTESIAN COORDINATES \(BOHR\)/;
    }
    $iatom = 0;
    while (<STDIN>)
    {
      chomp;
      (undef, $charge, $x, $y, $z) = split;
      last if ($charge !~ /\d\.0/);
      ($x, $y, $z) = map {$_*$BORHTOA} ($x, $y, $z);
      $atom[$iatom++] = join (' ', $charge, $x, $y, $z);
      $label{$charge} = 'a';
      $natoms{$charge}++ unless $natoms;
    }
    last RPATH unless $iatom;      # no more coordinates
    $natoms = $iatom unless $natoms;
    if ($opt_o eq 'i')
    {
      @atom = sort bycharge @atom;
      $ilabel = 0;
      foreach $charge (sort bycharge keys %label)
      {
        $label{$charge} .= ++$ilabel;
      }
      $nlabels = $ilabel;
    }
    if ($rpath > 0)        # evtl create output file
    {
      if ($ARGV[1])
      {
        $point = sprintf "%03d", $rpath;
        open (STDOUT, ">$ARGV[1]$point") ||
          die "Cannot open output file $ARGV[1]$point";
      }
      geniglo    () if $opt_o eq 'i';    # generate IGLO output file
      gencloppa  () if $opt_o eq 'c';    # generate CLOPPA output file
      gendalton  () if $opt_o eq 'd';    # generate DALTON output file
      genpriroda () if $opt_o eq 'p';    # generate PRIRODA output file
      genmolden  () if $opt_o eq 'm';    # generate MOLDEN xyz file
    }
  }
}

# Subroutines to generate output files...

# generate IGLO input file
sub geniglo
{
  my ($iatom, $charge, $x, $y, $z, $orb, $scfparam, @basis);
  
  # print title
  if ($rpath < 2)        # concatenate title only once
  {
    $title .= ', ' if $title;
    $title .= $gsymm;
    $title .= ', ' if $title;
    $title .= $basisname;
    $title .= ', ' if $title;
    $title .= 'IGLO';
  }
  print " ***, $title\n";

  # print memory
  printf "\n %-30s! memory requirement\n", $memory;

  # print usual commands
  print "\n int,nosort\n";

  # print symmetry
  printf " %-30s! %s symmetry\n\n", $symm, $gsymm;

  # print geometry
  for ($iatom=0; $iatom<$natoms; $iatom++)
  {
    ($charge, $x, $y, $z) = split (' ', $atom[$iatom]);
    printf " %s, %-2s, %10.6f, %10.6f, %10.6f, ang\n",
      $label{$charge}, $atname[$charge], $x, $y, $z;
  }
  print "\n";

  # print fallback ECP if any
  foreach $charge (sort bycharge keys %label)
  {
    if (! $basis{$atname[$charge]} && $basisecp{$atname[$charge]})
    {
      printf " ecp,%s,%s ! only ECP basis for %s\n",
        substr($label{$charge},1), $basisecp{$atname[$charge]},
        $atname[$charge];
      printf STDERR "!!! Only ECP basis for %s !!!\n", $atname[$charge];
    }
  }

  # print basis or fallback basis what appropriate
  foreach $charge (sort bycharge keys %label)
  {
    if ($basis{$atname[$charge]})
    {
      @basis = split (' ', $basis{$atname[$charge]});
    }
    else
    {
      @basis = split (' ', $basisfb{$atname[$charge]});
    }
    printf " ";
    for ($orb=0; $orb<$#basis; $orb++)
    {
      printf "%s,%s,%s;",
        $orbname[$orb], substr($label{$charge},1), $basis[$orb];
    }
    printf "%s,%s", $orbname[$orb], substr($label{$charge},1);
    if ($basis{$atname[$charge]})
    {
      printf ",%s\n", $basis[$orb];
    }
    else
    {
      if ($basis[$orb])
      {
  printf ",%s", $basis[$orb];
      }
      printf " ! no %s for %s\n", $basisname, $atname[$charge];
      printf STDERR "!!! No %s for %s !!!\n", $basisname, $atname[$charge];
    }
  }

  # print usual commands
  print " -\n";
  printf " %-30s! forget about small integrals\n", "thr,9,8";
  print "\n iglo\n";

  #print charge
  $icharg = sprintf "charge,%d", $icharg;
  printf " %-30s! charge of molecule\n", $icharg;

  # print usual commands
  $scfparam = sprintf "scfparam,50,4,3,%d,7", $istart;
  printf " %-30s! if not convergent, try scfparam,50,4,-1,%d,7\n",
    $scfparam, !$istart;
}

# generate CLOPPA input file
sub gencloppa
{
  my ($iatom, $charge, $x, $y, $z);

  # print comment
  print "CLOPPA-INDO\n";

  # print directives
  print "INDO XYZ NMRJ FERM SPIN ORBI END\n";

  # print options
  if ($icharg)
  {
    $icharg = sprintf ("CHARGE =%d ", $icharg) if $icharg;
  }
  else
  {
    $icharg = '';
  }
  printf "%sDENMATR=0 BINDING=0 BONDDEN=0 VECTORS=0 END\n", $icharg;

  # print title
  print "$title\n";

  # print geometry
  for ($iatom=0; $iatom<$natoms; $iatom++)
  {
    ($charge, $x, $y, $z) = split (' ', $atom[$iatom]);
    printf "%3d %10.6f %10.6f %10.6f\n", $charge, $x, $y, $z;
  }
  printf "%3d %10.6f %10.6f %10.6f\n", 0, 0, 0, 0;

  # print atom numbers
  printf "%d", $natoms;
  for ($iatom=0; $iatom<$natoms; $iatom++)
  {
    printf " %d", $iatom+1;
  }
  print "\n";
}

# generate DALTON MOL input file
sub gendalton
{
  my (@gensymm, $iatom, $charge, $n, $x, $y, $z, $name);
  
  print "ATOMBASIS\n";
  
  # print title
  print "$title\n";
  print "$gsymm symmetry, basis $opt_b\n";

  # print control line
  @gensymm = split (',', $symm);
  printf " %4d%3d%2d%3s%3s%3s1\n",
    scalar keys %label, $icharg, scalar @gensymm,
    "\U$gensymm[0]", "\U$gensymm[1]", "\U$gensymm[2]";

  # print geometry
  foreach $charge (sort bycharge keys %label)
  {
    printf "      %3d.%5d %s\n", $charge, $natoms{$charge}, $opt_b;
    for ($iatom=0; $iatom<$natoms; $iatom++)
    {
      ($n, $x, $y, $z, $name) = split (' ', $atom[$iatom]);
      next unless $n == $charge;
      $name = $atname[$charge].($iatom+1) unless $opt_s;
      $name = $atname[$charge].($iatom+1) unless $name;
      printf "%-4.4s %10.6f %10.6f %10.6f\n", $name, $x, $y, $z;
    }
  }
}

# generate PRIRODA input file
sub genpriroda
{
  my ($iatom, $charge, $x, $y, $z);

  # print comment
  print "Input for Priroda from $ARGV[0] by g2i";
  print "     $title" if ($title);
  # print usual commands
  print '
$system mem=64 disk=1 $end
$control
 task=nmr',"\n";
  # print basis
  if ($opt_b) {
    print "basis=$opt_b.bas";
  }
  else {
    print "basis=3z.bas";
  }
  # print usual commands
  print '
 print=+charges
$end
$grid accur=1e-8 $end
$optimize
 steps=100
 tol=1e-5
$end
$molecule',"\n ";
  # print charge
  print "charge=$icharg " if ($icharg); 
  print "mult=1
 cartesian\n";
  # print geometry
  for ($iatom=0; $iatom<$natoms; $iatom++)
  {
    ($charge, $x, $y, $z) = split (' ', $atom[$iatom]);
    printf "%3d %12.8f %12.8f %12.8f\n", $charge, $x, $y, $z;
  }
  print '$end',"\n";
}

sub genmolden
{
  my ($iatom, $charge, $x, $y, $z);

  # 1-st line
  print " $natoms\n";
  # 2-nd line
  print " Energy $energy" if ($energy);
  print "   POINT $rpath" if ($rpath>1);
  print "   $title" if ($title);
  print "\n";
  # print geometry
  for ($iatom=0; $iatom<$natoms; $iatom++) {
    ($charge, $x, $y, $z) = split (' ', $atom[$iatom]);
    printf " %-2s %12.8f %12.8f %12.8f\n", $atname[$charge], $x, $y, $z;
  }
}

# Utility subroutines...

# reverse sort by atom charge
sub bycharge
{
  my ($chargea) = split (' ', $a);
  my ($chargeb) = split (' ', $b);
  $chargeb <=> $chargea;
}

# Actual program start...

# get options
getopts ('b:hm:o:s');

# short help
if ($opt_h)
{
  print '
Извлекает из GAMESS and Priroda output files последнюю геометрию
и печатает ее в MOLDEN-формате, либо, если заданы соответствующие опции,
Converts GAMESS and Priroda output files
into IGLO, CLOPPA, Priroda and Dalton input files

Usage: g2i [-b <basis>] [-m <memory>] [-o <type>] [-s] [<infile>] [<outfile>]

<type> of output format is MOLDEN by default

If <type> is C (case insensitive) - CLOPPA
   <basis>, <memory> and -s are irrelevant

If <type> is D (case insensitive) - DALTON
   <basis> is the name of the basis description file
   <memory> is irrelevant
   If -s is given, symbolic atom names from <infile> are used

If <type> is I (case insensitive) - IGLO
   <basis> may be DZ (default), I, II, III, IV (case insensitive) or 1, 2, 3, 4
   <memory> may be given in k, K, m, M or 8-byte words
   -s is irrelevant

If <type> is P (case insensitive) - Priroda
   <basis> is the name of the basis description file (3z is default)
   <memory> and -s are irrelevant

If <type> is M (case insensitive) - xyz in MOLDEN format
   <basis>, <memory> and -s are irrelevant

<infile> and <outfile> are STDIN and STDOUT by default

The program knows GAMESS and Priroda file formats and recognizes them on 
input automatically. Input may be multipoint (SCAN or IRC) or concatenated 
GAMESS/Priroda files. In these cases output includes every point, and 
<outfile> is incremented',"\n";
  exit 0;
}

# type of output
$opt_o = "\L$opt_o";
$opt_o = 'm' if (! $opt_o);

# basis name
if (! $opt_b)
{
  if ($opt_o eq 'i')
  {
    $opt_b = 'DZ';
  }
  elsif ($opt_o eq 'd')
  {
    $opt_b = '6-31G';
  }
}
if (! $opt_b || "\L$opt_b" eq 'dz')
{
  %basis = %basisdz;
  $basisname = 'basis DZ';
  $istart = 1;
}
elsif ($opt_b eq 1 || "\L$opt_b" eq 'i')
{
  %basis = %basis1;
  $basisname = 'basis I';
  $istart = 1;
}
elsif ($opt_b eq 2 || "\L$opt_b" eq 'ii')
{
  %basis = %basis2;
  $basisname = 'basis II';
  $istart = 1;
}
elsif ($opt_b eq 3 || "\L$opt_b" eq 'iii')
{
  %basis = %basis3;
  $basisname = 'basis III';
  $istart = 0;
}
elsif ($opt_b eq 4 || "\L$opt_b" eq 'iv')
{
  %basis = %basis4;
  $basisname = 'basis IV';
  $istart = 0;
}
elsif ($opt_o eq 'i')
{
  printf STDERR "!!! Basis %s not recognized !!!\n", $opt_b;
  %basis = %basisdz;
  $basisname = 'basis DZ';
  $istart = 1;
}

# scan input file
if ($ARGV[0])
{
  open (STDIN, $ARGV[0]) || die "Cannot open input file $ARGV[0]";
}

# determine the creator of the file
$outtype = 'unknown';
TYPE:while (<STDIN>)
{
  chomp;
  if (/^\s*\*{30,}$/)
  {
    do
    {
      $_ = <STDIN>;
      chomp;
      last unless /^\s*\*\s+GAMESS\s+VERSION\s+=\s+.+\s+\*$/;
      $_ = <STDIN>;
      chomp;
      last unless /^\s*\*\s+FROM\s+IOWA\s+STATE\s+UNIVERSITY\s+\*$/;
      $outtype = 'gamess';
      last TYPE;
    }
  }
  if (/^\s*program\s+Priroda\s+version\s+.+$/)
  {
    do
    {
      $_ = <STDIN>;
      chomp;
      last unless /^\s*copyright\s+\(c\)\s+.+\s+Dimitri\s+Laikov$/;
      $outtype = 'priroda';
      last TYPE;
    }
  }
  if (/^\*\*\*\*\* BEGIN IRC INFORMATION PACKET \*\*\*\*\*$/)
  {
    do
    {
      $outtype = 'gamess_scr_irc';
      seek (STDIN, 0, 0);
      last TYPE;
    }
  }
}
if ($outtype eq 'unknown')
{
  seek (STDIN, 0, 0);  # no stamp found - rewind and set default to gamess
  $outtype = 'gamess';
}

# interpret input file
$rpath = 0;          # no reaction path found yet
readgamess  () if $outtype eq 'gamess';    # try read GAMESS input file
readpriroda () if $outtype eq 'priroda';  # try read Priroda input file
readscrirc  () if $outtype eq 'gamess_scr_irc';  # try read gamess_scr_irc input file

# create output file if not yet done for stationary points
if ($ARGV[1])
{
  open (STDOUT, ">$ARGV[1]") || die "Cannot open output file $ARGV[1]";
}
geniglo    () if $opt_o eq 'i';      # generate IGLO output file
gencloppa  () if $opt_o eq 'c';      # generate CLOPPA output file
gendalton  () if $opt_o eq 'd';      # generate DALTON output file
genpriroda () if $opt_o eq 'p';      # generate Priroda output file
genmolden  () if $opt_o eq 'm';      # generate MOLDEN xyz file

# End of file.