#!/usr/bin/perl -ws

our ($h,$help);
(my $program = $0) =~ s/^.*[\/\\]//;

if ($h || $help) {
  print <<H;

Gaussian fchk to xyz convertor (for IRC tasks)

Usage: $0 checkpoint.fchk > IRC.xyz

H
  exit;
}

die "Usage: $0 checkpoint.fchk > IRC.xyz\n" unless @ARGV==1;

# 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
       );

#use Data::Dump 'pp';

#my %chk = read_fchk('STDIN', 'Number of atoms', 'Atomic numbers');
my %chk = read_fchk($ARGV[0]);
#pp keys(%chk);

my $N = $chk{'Number of atoms'};
my @atom = map {$atname[$_]} @{$chk{'Atomic numbers'}};
my $IRC_steps = $chk{'IRC Number of geometries'}[0];

my $charge = $chk{Charge} if $chk{Charge} != 0;
my $mult = $chk{Multiplicity} if $chk{Multiplicity} != 1;

my @energy = @{$chk{'IRC point       1 Results for each geome'}};
my @coord = @{$chk{'IRC point       1 Geometries'}};
my @grad = @{$chk{'IRC point       1 Gradient at each geome'}};

my @mols;
if (grep {m/^IRC/} keys %chk) {
  for (my $i=0; $i<$IRC_steps; $i++) {
    $mols[$i][0]{Energy} = $energy[2*$i];
    $mols[$i][0]{Rx_Coord} = $energy[2*$i+1];
    $mols[$i][0]{Charge} = $charge if $charge;
    $mols[$i][0]{Mult} = $mult if $mult;
    $mols[$i][0]{Gradient} = sqrt(sum(map {$_*$_} @grad[$i*$N*3..($i+1)*$N*3-1]));
    for (my $j=1; $j<=$N; $j++) {
      $mols[$i][$j][0] = $atom[$j-1];
      $mols[$i][$j][1] = $coord[3*($N*$i+$j-1)]*0.529177249;
      $mols[$i][$j][2] = $coord[3*($N*$i+$j-1)+1]*0.529177249;
      $mols[$i][$j][3] = $coord[3*($N*$i+$j-1)+2]*0.529177249;
    }
  }

  my (@mol_p,@mol_m);
  foreach my $mol (@mols) {
    if ($mol->[0]{Rx_Coord}<0) {
      push @mol_m, $mol;
    } 
    else {
      push @mol_p, $mol;
    }
  }
  @mols = (reverse(@mol_m), @mol_p);
}
else {
  warn "No IRC data found\n";
}

#pp @mols;
write_molden(@mols) if @mols;


sub read_fchk {
  my %chk;
  my $file = shift;
  open L, '<', $file or die "Can't open $file: $!\n";
  while (<L>) {
    my ($item,$IR,$is_array,$rest) = unpack("A43 A1 x3 A2 A*", $_);
    next if $item =~ /^ /;
    next if $IR !~ /[IR]/;
    $item =~ s/\s+$//;
    if ($is_array eq 'N=') {
      my $n = 0+$rest;
      my @a;
      while (@a < $n) {
        push @a, split(' ', <L>);
      }
      die "Bad for $item\n" if @a != $n;
      $chk{$item} = \@a;
    }
    else {
      $rest =~ s/^\s+//;
      $chk{$item} = $rest;
    }
  }
  close L;
  return %chk;
}

sub write_molden {
  my $fh;
  unless (ref $_[-1]) {
    my $file = pop @_;
    open $fh, '>', $file or do {warn "Can't write to $file: $!\n"; return};
  }
  else {
    $fh = 'STDOUT';
  }
  foreach my $mol (@_) {;
    my $N = $#{$mol};
    print $fh " $N\n";
    printf $fh " Energy %.8f", $mol->[0]{Energy} if defined $mol->[0]{Energy};
    printf $fh "  Rx_Coord %.4f", $mol->[0]{Rx_Coord} if defined $mol->[0]{Rx_Coord};
    printf $fh "  Gradient %.8f", $mol->[0]{Gradient} if defined $mol->[0]{Gradient};
    print  $fh "  Charge $mol->[0]{Charge}" if defined $mol->[0]{Charge};
    print  $fh "  Mult $mol->[0]{Mult}" if defined $mol->[0]{Mult};
    print  $fh "  Symmetry $mol->[0]{Symmetry}" if $mol->[0]{Symmetry};
    print  $fh "\n";
    for (my $i=1; $i<=$N; $i++) {
      printf $fh " %-2s %12.8f %12.8f %12.8f\n", @{$mol->[$i]};
    }
  }
  close $fh if $fh ne 'STDOUT';
}

sub sum {
  my $sum = 0;
  foreach (@_) {
    $sum += $_;
  }
  return $sum;
}