#!/usr/bin/perl -ws

#use Data::Dump qw(pp);
use vars qw($h);
(my $program = $0) =~ s/^.*[\/\\]//;

if ($h) {
  print <<HELP;

Usage: $program files.xyzppm
       cat files.xyzppm | $program
Зависимости: perl, symmetry (не обязательно)

Печать таблицы хим. сдвигов, подходящей для импорта в для Excel, из *.xyzppm

Автоматически будут усреднены и объединены протоны метилов.

HELP
  exit;
}

#$symmetry = '/usr/local/bin/symmetry';
#unless (-x $symmetry) {
#  warn "Program symmetry not found\n";
#  undef $symmetry;
#}

# Вгоняем таблицу Менделеева
@ATOM = 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
           );
for ($i=0; $i<@ATOM; $i++) { $ATOM{$ATOM[$i]} = $i };

my $prev_head = '';

while (&read_molden) {
  my @new_atom = map {$atom[$_].$_} 1..$#atom;
  unshift @new_atom, undef;
  #pp @new_atom;
  
  # Усредняем протоны метильных групп
  for ($i=1; $i<=$N; $i++) {
    if ( $atom[$i] eq 'C' ) {
      $CH[$i] = [];
      for ($j=1; $j<=$N; $j++) {
        if ( $atom[$j] eq 'H' &&
             ($x[$i]-$x[$j])**2+($y[$i]-$y[$j])**2+($z[$i]-$z[$j])**2 < 1.44
            ) {
          push @{$CH[$i]}, $j;
        }
      }
    }
  }

  for ($i=1; $i<=$N; $i++) {
    if ($atom[$i] eq 'C' && @{$CH[$i]}==3) {
      my $aver = ($ppm[$CH[$i][0]]+$ppm[$CH[$i][1]]+$ppm[$CH[$i][2]])/3;
      $ppm[$CH[$i][0]] = $ppm[$CH[$i][1]] = $ppm[$CH[$i][2]] = $aver;
      $new_atom[$CH[$i][0]] .= ",$new_atom[$CH[$i][1]],$new_atom[$CH[$i][2]]";
      $new_atom[$CH[$i][1]] = $new_atom[$CH[$i][2]] = '';
    }
  }
  #pp @CH;

  my $head = "file\tEnergy";
  my $str = "$ARGV\t$energy";
  for ($i=1; $i<@new_atom; $i++) {
    next unless $new_atom[$i];
    next unless $ppm[$i];
    $head .= "\t$new_atom[$i]";
    $str .= sprintf "\t%.3f", $ppm[$i];
  }
  #pp @new_atom;
  print "$head\n" if $head ne $prev_head;
  print "$str\n";
  $prev_head = $head;

#  if ($symmetry) {
#    open TINP, ">symm.TINP" || warn "Can't write temporary file: $!\n" , next;
#    print TINP " $N\n";
#    for ($i=1; $i<=$N; $i++) {
#      print TINP " $ATOM{$atom[$i]}   $x[$i]  $y[$i]  $z[$i]\n";
#    }
#    close TINP;
#
#    $generator = undef;
#    open OUT, "$symmetry -final 0.0009999 symm.TINP |";
#    while (<OUT>) {
#      if (/^Molecule has the following abelian symmetry generators: (\S+)/) {
#        $generator = $1;
#        last;
#      }
#    }
#    while (<OUT>) {
#      if (m/^Abelian symmetrized coordinates of all atoms$/) {
#        <OUT>;
#        $N = <OUT>; chomp $N;
#        for ($i=1; $i<=$N; $i++) {
#          ($natom,$x[$i],$y[$i],$z[$i]) = split ' ', <OUT>;
#          $atom[$i] = $ATOM[$natom];
#        }
#      }
#    }
#    close OUT;
#    unlink 'symm.TINP';
#
#    @symm_atoms = symm_atoms($generator);
#    pp @symm_atoms;
#    SYMM_ATOMS:
#    foreach (@symm_atoms) {
#      #print "@{$_}\n";
#      my $sum = 0;
#      foreach (@{$_}) {
#        next SYMM_ATOMS if $ppm[$_] eq '';
#        $sum += $ppm[$_];
#      }
#      my $aver = $sum/@{$_};
#      $new_atom[$_->[0]] = join ',', map {$new_atom[$_]} @$_;
#      $ppm[shift @{$_}] = $aver;
#      foreach (@{$_}) {
#        $ppm[$_] = '';
#        $new_atom[$_] = '';
#      }
#    
#    }
#  }
  
  #&write_molden;
}

#sub symm_atoms {
#  my $generator = shift;
#  my @symm_atoms;
#  unless ($generator) {
#    @symm_atoms = map {[$_]} 1..$N;
#    return @symm_atoms;
#  }
#  my %condition = (
#    'z'        => '$x[$i]==$x[$j] && $y[$i]==$y[$j] && $z[$i]==-$z[$j]',
#    'x,y'      => 'abs($x[$i])==abs($x[$j]) && 
#                   abs($y[$i])==abs($y[$j]) && 
#                   $z[$i]==$z[$j]',
#    'xy'       => '$x[$i]==-$x[$j] && $y[$i]==-$y[$j] && $z[$i]==$z[$j]',
#    'x,y,z'    => 'abs($x[$i])==abs($x[$j]) && 
#                   abs($y[$i])==abs($y[$j]) && 
#                   abs($z[$i])==abs($z[$j])',
#    'xyz'      => '$x[$i]==-$x[$j] && $y[$i]==-$y[$j] && $z[$i]==-$z[$j]',
#    'xy,z'     => '($x[$i]==-$x[$j] && $y[$i]==-$y[$j] && $z[$i]==$z[$j]) || 
#                   ($x[$i]==$x[$j] && $y[$i]==$y[$j] && $z[$i]==-$z[$j])',
#    'xy,xz,yz' => '($x[$i]==-$x[$j] && $y[$i]==-$y[$j] && $z[$i]==$z[$j]) ||
#                   ($x[$i]==-$x[$j] && $y[$i]==$y[$j] && $z[$i]==-$z[$j]) ||
#                   ($x[$i]==$x[$j] && $y[$i]==-$y[$j] && $z[$i]==-$z[$j])'
#  );
#  my @at = @atom;
#  for ($i=1; $i<=$N; $i++) {
#    next unless $at[$i];
#    push @symm_atoms, [($i)];
#    #$at[$i] = undef;
#    for ($j=$i+1; $j<=$N; $j++) {
#      next unless $at[$j];
#      if (exists $condition{$generator} && eval "$condition{$generator}") {
#          $at[$j] = undef;
#          push @{$symm_atoms[-1]}, $j;
#      }
#    }
#  }
#  return @symm_atoms;
#}

sub read_molden {
  return undef if eof();
  # 1-st string
  if (<>=~/^\s*(\d+)\s*$/) {
    ($N,$energy,@atom,@x,@y,@z,@ppm) = undef;
    $N = $1;
  } else {
    return undef;
  }
  # 2-nd string
  ($energy) = <>=~/(-?\d+\.\d+)/;
  # Geometry
  for ($i=1;$i<=$N;$i++) {
    ($atom[$i],$x[$i],$y[$i],$z[$i],$ppm[$i],undef) = split ' ', <>;
    return undef unless $atom[$i]=~/^\w\w?/ &&
                        $x[$i]=~/^-?\d+\.\d+/ &&
                        $y[$i]=~/^-?\d+\.\d+/ &&
                        $z[$i]=~/^-?\d+\.\d+/;
  }
  return $N;
}

sub write_molden {
  if (@atom && @x && @y && @z) {
    print " $#atom\n";
    print " Energy $energy" if $energy;
    print "\n";
    for ($i=1; $i<=$#atom; $i++) {
      if ($atom[$i] && $x[$i] && $y[$i] && $z[$i]) {
        printf " %-2s %12.8f %12.8f %12.8f",$atom[$i],$x[$i],$y[$i],$z[$i];
        printf $atom[$i] eq 'H' ? " %10.3f" : " %9.2f"  , $ppm[$i] if $ppm[$i];
        print "\n";
      }
    }
  }
}