#!/usr/bin/perl -ws
#use Data::Dump;

our ($h,$help,$join);
if ($h || $help) {
  (my $program = $0) =~ s/^.*[\/\\]//;
  print <<HELP;
Usage: $program file.xyz > dihedrals
Generate list of dihedrals 
(dihedral per line, numbers of atoms through spaces)
HELP
  exit;
}


# Ковалентные радиусы из Dalton Trans., 2008, 2832-2838
our %radius = qw(
 H 0.31  He 0.28  Li 1.28  Be 0.96   B 0.84   C 0.73   N 0.71   O 0.66   F 0.57  
Ne 0.58  Na 1.66  Mg 1.41  Al 1.21  Si 1.11   P 1.07   S 1.05  Cl 1.02  Ar 1.06  
 K 2.03  Ca 1.76  Sc 1.70  Ti 1.60   V 1.53  Cr 1.39  Mn 1.50  Fe 1.44  Co 1.38   
Ni 1.24  Cu 1.32  Zn 1.22  Ga 1.22  Ge 1.20  As 1.19  Se 1.20  Br 1.20  Kr 1.16  
Rb 2.20  Sr 1.95   Y 1.90  Zr 1.75  Nb 1.64  Mo 1.54  Tc 1.47  Ru 1.46  Rh 1.42   
Pd 1.39  Ag 1.45  Cd 1.44  In 1.42  Sn 1.39  Sb 1.39  Te 1.38   I 1.39  Xe 1.40   
Cs 2.44  Ba 2.15  La 2.07  Ce 2.04  Pr 2.03  Nd 2.01  Pm 1.99  Sm 1.98  Eu 1.98   
Gd 1.96  Tb 1.94  Dy 1.92  Ho 1.92  Er 1.89  Tm 1.90  Yb 1.87  Lu 1.87  Hf 1.75   
Ta 1.70   W 1.62  Re 1.51  Os 1.44  Ir 1.41  Pt 1.36  Au 1.36  Hg 1.32  Tl 1.45  
Pb 1.46  Bi 1.48  Po 1.40  At 1.50  Rn 1.50  Fr 2.60  Ra 2.21  Ac 2.15  Th 2.06  
Pa 2.00   U 1.96  Np 1.90  Pu 1.87  Am 1.80  Cm 1.69  
);

my $file = $ARGV[0];
my ($mol) = read_molden($file);
my @dihs = get_dihs_list($mol);
print "$_\n" foreach (@dihs);

sub get_dihs_list {
  my $mol = shift;
  my $N = @$mol;
  my @bonds; # i-j
  my @valence; # of i and j
  my $f = 1.3;
  for (my $i=1; $i<$N; $i++) {
    my ($ai,$xi,$yi,$zi) = @{$mol->[$i]};
    die "No covalent radius for $ai\n" unless exists $radius{$ai};
    for (my $j=$i+1; $j<$N; $j++) {
      my ($aj,$xj,$yj,$zj) = @{$mol->[$j]};
      my $d = sqrt(($xi-$xj)**2 + ($yi-$yj)**2 + ($zi-$zj)**2);
      if ($d < ($radius{$ai}+$radius{$aj})*$f ) {
        push @bonds, [$i,$j];
        push @{$valence[$i]}, $j;
        push @{$valence[$j]}, $i;
      }
    }
  }
  
  my @dihs_list; # k-i-j-l
  foreach my $bond (@bonds) {
    my ($i,$j) = @$bond;
    my @k = @{$valence[$i]};
    my @l = @{$valence[$j]};
    
    # Exclude bonds with monovalent i or j
    next if @k==1 or @l==1;
    
    # Delete j from k and i from l
    @k = grep {$_ != $j} @k;
    @l = grep {$_ != $i} @l;
    
    # Exclude H-C-j-l and k-i-C-H
    if ($mol->[$i][0] eq 'C') {
      @k = grep {$mol->[$_][0] ne 'H'} @k;
    }
    if ($mol->[$j][0] eq 'C') {
      @l = grep {$mol->[$_][0] ne 'H'} @l;
    }
    next if ! @k or ! @l; # Terminal Me
    
    #print "[@k]-$i-$j-[@l]\n";
    # Expand @k and @l
    foreach my $k (@k) {
      foreach my $l (@l) {
        push @dihs_list, "$k $i $j $l";
      }
    }
  }
  return @dihs_list unless $join;
  
  LOOP:
  while (1) {
    #print join "\n", @dihs_list, "\n";
    foreach my $D (@dihs_list) {
      next unless $D;
      (my $KIJ) = $D=~/^(\d+ \d+ \d+)/;
      (my $IJL) = $D=~/(\d+ \d+ \d+)$/;
      foreach my $d (@dihs_list) {
        next unless $d;
        next if $d eq $D;
        (my $kij) = $d=~/^(\d+ \d+ \d+)/;
        (my $ijl) = $d=~/(\d+ \d+ \d+)$/;
        if ($kij eq $IJL) {
          #print "  $D,$d:\n";
          $d =~ s/$kij//;
          $D = $D . $d;
          $d = '';
          next LOOP;
        }
        if ($ijl eq $KIJ) {
          #print "  $D,$d:\n";
          $d =~ s/$ijl//;
          $D = $d . $D;
          $d = '';
          next LOOP;
        }
        $d = join ' ', reverse split ' ', $d;
        next if $d eq $D;
        ($kij) = $d=~/^(\d+ \d+ \d+)/;
        ($ijl) = $d=~/(\d+ \d+ \d+)$/;
        if ($kij eq $IJL) {
          #print "  $D,$d:\n";
          $d =~ s/$kij//;
          $D = $D . $d;
          $d = '';
          next LOOP;
        }
        if ($ijl eq $KIJ) {
          #print "  $D,$d:\n";
          $d =~ s/$ijl//;
          $D = $d . $D;
          $d = '';
          next LOOP;
        }
      }
    }
    last LOOP;
  }
  return grep {$_} @dihs_list;
}


sub read_molden {
 ############################################################################
 ## Молекула -- ссылка на массив, 0-й элемент -- свойства, следующие - атомы.
 ## Свойства -- ссылка на хэш с ключами Energy, Symmetry
 ## Атом -- ссылка на массив [atom, x, y, z, ppm] (ppm может не быть).
 ##
 ## Читает xyz. Параметры - имена xyz-файлов. Если параметров нет, то <>.
 ## Возвращает массив найденных молекул.
 ############################################################################
  local @ARGV = @_ ? @_ : @ARGV;
  my $num = qr/-?\d+(?:\.\d+)?/;
  my @mols;
  my $line;
  LOOP:
  while ($line || defined($line = <>)) {
    #print $line;
    if ($line =~/^\s*(\d+)\s*$/) {
      my @mol;
      my $N = $1;
      last LOOP if eof();
      next LOOP if eof(ARGV);
      $line = <>;
      ($mol[0]{'Energy'}) = $line =~ /(?:^|)\s($num)(?:\s|$)/o;
      ($mol[0]{'Symmetry'}) = $line =~ /symm\S*\s+(\S+)/i;
      for (my $i=1; $i<=$N; $i++) {
        last LOOP if eof();
        next LOOP if eof(ARGV);
        $line = <>;
        #print $line;
        if ($line =~ /^\s*([A-Z]{1,2})\s+($num)\s+($num)\s+($num)\s*(.*)/io) {
          $mol[$i] = [$1,$2,$3,$4,$5];
        } else {
          next LOOP;
        }
      }
      push @mols, \@mol;
      last LOOP if eof();
    } else {
      undef $line;
    }
  }
  return @mols;
}