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

our ($h,$help);
(my $program = $0) =~ s/^.*[\/\\]//;
if ($h || $help) {
  print "
  Compare spatial structures for each of files.xyz with each of FILES.xyz.
  Comparison is independent of order of atom's numeration.
  
  Usage: $program files.xyz vs. FILES.xyz
  'vs.' divides comparable molecules (from files.xyz) 
  and reference ones (from FILES.xyz)
  
  Another usage: cat files.xyz | $program FILES.xyz
  f.e. xyz -n=1,-1 IRC.xyz | $program *.O.xyz
  
  If file.xyz contains several molecules (concatenated xyz) 
  then each molecule is named as 'file_n'. Piped molecules is named as 'pipe_n'.
  
  -chiral  optical isomers are not identical
  -RMS=0.1  molecules are not identical if RMSD > 0.1
  -W  all atom weight to be equal 1 (default RMS is in sqrt(amu)*anstrom units).
      Don't take into account atom types for -dist.
  -H  don't take into account hydrogens
  -all  compare all files with each other: $program -all FILES.xyz
  -pair compare 1st from files.xyz with 1st from FILES.xyz, 2nd with 2nd, etc
  -dih='i,j,k,l i1,j1,k1,l1, ...' compare by dihedrals, default -RMS=30 deg. 
  -dih=i,j,k,l,m,n,i,j,k  shorthand for sequential dihedrals (e.g. for cycles)
  -dihf=dihedrals  compare by dihedrals provided in 'dihedrals' file
     (dihedral per line, numbers of atoms through spaces)
     Automatic genaration of dihedrals: dih_list file.xyz > dihedrals
  -dist  use modified Grigoryan-Springborg algorithm (distance matrix) for 
         comparison. -chiral doesn’t affect
  -renum  $program -renum *.xyz  OR  cat *.xyz | $program -renum
          Renumbers the structures (each by the previous one) and 
          prints xyz to stdout.
  -verbose  more detailed output
";
  exit;
}

our ($chiral,$RMS,$verbose,$all,$pair,$dih,$dihf,$H,$W,$dist,$renum);

if (defined $dihf) {
  $dihf = 'dihedrals' if $dihf eq '1';
  die "No dihedrals $dihf\n" unless -f $dihf;
  open L, '<', $dihf or die "Can't open $dihf: $!\n";
  while (<L>) {
    chomp;
    s/^\D+//;
    s/\D+$//;
    $dih .= ' ' . join(',', split /\D+/, $_);
  }
  close L;
}

# Массы изотопов
our %massa = qw(XX 0.000
H  1.007947    He 4.0026022   Li 6.9412      Be 9.0121823   B  10.8117
C  12.01078    N  14.00672    O  15.99943    F  18.99840325 Ne 20.17976
Na 22.9897702  Mg 24.30506    Al 26.9815382  Si 28.08553    P  30.9737612
S  32.0655     Cl 35.4532     Ar 39.9481     K  39.09831    Ca 40.0784
Sc 44.9559108  Ti 47.8671     V  50.94151    Cr 51.99616    Mn 54.9380499
Fe 55.8452     Co 58.9332009  Ni 58.69342    Cu 63.5463     Zn 65.4094
Ga 69.7231     Ge 72.641      As 74.921602   Se 78.963      Br 79.9041
Kr 83.7982     Rb 85.46783    Sr 87.621      Y  88.905852   Zr 91.2242
Nb 92.906382   Mo 95.942      Ru 101.072     Rh 102.905502  Pd 106.421
Ag 107.86822   Cd 112.4118    In 114.8183    Sn 118.7107    Sb 121.7601
Te 127.603     I  126.904473  Xe 131.2936    Cs 132.905452  Ba 137.3277
La 138.90552   Ce 140.1161    Pr 140.907652  Nd 144.243     Sm 150.363
Eu 151.9641    Gd 157.253     Tb 158.925342  Dy 162.5001    Ho 164.930322
Er 167.2593    Tm 168.934212  Yb 173.043     Lu 174.9671    Hf 178.492
Ta 180.94791   W  183.841     Re 186.2071    Os 190.233     Ir 192.2173
Pt 195.0782    Au 196.966552  Hg 200.592     Tl 204.38332   Pb 207.21
Bi 208.980382  Th 232.03811   Pa 231.035882  U  238.028913
);
# Ковалентные радиусы из 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  
);

local $| = 1;

# Periodic table
our @element = 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 Rf Db Sg Bh Hs Mt Ds Rg
              );
our %element;
@element{@element} = 0..$#element;


if ($dih && $dist) {
  die "Inconsistent options -dih and -dist\n";
}

$RMS ||= 0.1;
if ($dih) {
  $RMS ||= 30;
}
elsif ($dist) {
  $RMS /= 2;
}

if ($all && $pair) {
  die "Inconsistent options -all and -pair\n";
}

if ($W) {
  foreach (keys %massa) {
    $massa{$_} = 1;
  }
}

if ($verbose) {
  print "Chirality is taken into account.\n" if $chiral and !$dist;
  print "Chirality isn't taken into account.\n" unless $chiral;
  print "Compare all molecules with each other\n" if $all;
  print "Compare in pairs\n" if $pair;
  if ($dih) {
    print "Threshold max dihedral difference $RMS (degrees)\n";
  }
  elsif ($dist) {
    print "Threshold Grigoryan-Springborg criterion $RMS\n";
  }
  else {
    print $H ? "Threshold RMSD $RMS Angsrom without hydrogens\n"
             : "Threshold RMSD $RMS Angsrom*(amu)^1/2\n";
  }
}

my (@mols1,@mols2,$vs);
if (-p STDIN) { # If pipe
  die "Usage: cat files.xyz | $program FILES.xyz OR cat files.xyz | $program -all\n" if (! @ARGV && ! $all && ! $renum);
  my @files2 = @ARGV;
  undef @ARGV;
  @mols1 = read_molden();
  
  foreach my $file (@files2) {
    (my $name = $file) =~ s/\.xyz$//;
    my @mols = read_molden($file);
    next unless @mols;
    if (@mols > 1) {
      for (my $i=0; $i<@mols; $i++) {
        $mols[$i][0]{Name} = $name.'_'.($i+1);
        push(@mols2, $mols[$i]);
      }
    }
    else {
      $mols[0][0]{Name} = $name;
      push(@mols2, $mols[0]);
    }
  }
  
  my $count = 1;
  foreach my $mol (@mols1) {
    $mol->[0]{Name} = 'pipe_' . $count++;
  }
  $mols1[0][0]{Name} = 'pipe' if @mols1 == 1;
  #dd @mols1; print "\n\n\n"; dd @mols2; exit;
}
else { # No pipe
  if (! $all && ! $renum) {
    die "Usage: $program files.xyz vs. FILES.xyz OR $program -all files.xyz\n" unless grep {m/^vs\.$/} @ARGV;
  }
  foreach my $file (@ARGV) {
    if ($file =~ m/^vs\.$/) {
      $vs = 1;
      next;
    }
    (my $name = $file) =~ s/\.xyz$//;
    my @mols = read_molden($file);
    next unless @mols;
    if (@mols > 1) {
      for (my $i=0; $i<@mols; $i++) {
        $mols[$i][0]{Name} = $name.'_'.($i+1);
        $vs ? push(@mols2, $mols[$i]) : push(@mols1, $mols[$i]);
      }
    }
    else {
      $mols[0][0]{Name} = $name;
      $vs ? push(@mols2, $mols[0]) : push(@mols1, $mols[0]);
    }
  }
}

if ($all or $renum) {
  @mols2 = (@mols1,@mols2);
  die "No molecules for comparison\n" unless @mols2;
}
else {
  die "No molecules in comparable files\n" unless @mols1;
  die "No molecules in reference FILES\n" unless @mols2;
}
if ($pair) {
  die "Different numbers of molecules in files and FILES\n" if @mols1!=@mols2;
}

my @names = map {$_->[0]{Name}} @mols2;
delete_eq_end(\@names) if ! $pair;
#print "@names\n"; exit;
for (my $i=0; $i<@mols2; $i++) {
  $mols2[$i][0]{Name} = $names[$i];
}


#dd @mols1; print "\n\n"; dd @mols2;

if ($H) {
  # Atom numbers may be changed after H's removing, therefore
  die "-H and -dih are incompatible options\n" if $dih;
#  foreach my $mol (@mols1) {
#    @$mol = grep {! ref || ref ne 'ARRAY' || uc($_->[0]) ne 'H'} @$mol;
#  }
#  foreach my $mol (@mols2) {
#    @$mol = grep {! ref || ref ne 'ARRAY' || uc($_->[0]) ne 'H'} @$mol;
#  }
}

# Delete XX and H's
foreach my $m (@mols1,@mols2) {
  my @m1 = ($m->[0]);
  for (my $i=1; $i<@$m; $i++) {
    next if uc($m->[$i][0]) eq 'XX';
    next if $H && uc($m->[$i][0]) eq 'H';
    $m->[$i][0] = ucfirst(lc $m->[$i][0]);
    push @m1, $m->[$i];
  }
  $m = \@m1;
  # pp $m;
}

# Делаем проверку на число атомов. Нумеруем молекулы
if (! $pair) {
  for (my $i=0; $i<@mols2; $i++) {
    if (get_formula($mols2[$i]) ne get_formula($mols2[0])) {
      die "$mols2[$i][0]{Name} is not isomer of $mols1[0][0]{Name}. Die.\n";
    }
    $mols2[$i][0]{'Point'} = $i+1;
  }
  if (! $all) {
    for (my $i=0; $i<@mols1; $i++) {
      if (get_formula($mols1[$i]) ne get_formula($mols2[0])) {
        die "$mols1[$i][0]{Name} is not isomer of $mols1[0][0]{Name}. Die.\n";
      }
      $mols1[$i][0]{'Point'} = $i+1;
    }
  }
}

my $start = time();
if ($dist) {
  print "\nCalculation of distance matrix: " if $verbose;
  distance_matrix(@mols1,@mols2);
}
elsif ($dih) {
  print "\nCalculation of dihedrals: " if $verbose;
  calculate_dihedrals(@mols1,@mols2,$dih);
}
else {
  print "\nOrientation of molecules by inertia ellipsoid: " if $verbose;
  #my @eigen_num1 = tenzor(@mols1,@mols2);
  tenzor(@mols1,@mols2);
  #exit;
  #my @eigen_num2 = tenzor(@mols2);
}
print time()-$start, " s\n\n" if $verbose;

#exit;

print "Comparison:\n" if $verbose;
if ($all) {
  my (@nodes,@edges);
  for (my $i=0; $i<@mols2; $i++) {
    my $mol1 = $mols2[$i];
    push @nodes, $mol1->[0]{Name};
    for (my $j=$i+1; $j<@mols2; $j++) {
      my $mol2 = $mols2[$j];
      my $rms;
      if    ($dih)  { $rms = dih_rms($mol1,$mol2) }
      elsif ($dist) { $rms = GrigoryanSpringborg($mol1,$mol2) }
      else          { $rms = tenzor_sort($mol1,$mol2) }
      printf "$mol1->[0]{Name} vs. $mol2->[0]{Name}: diff=%.3f\n",$rms if $verbose;
      if ($rms < $RMS) {
        print "$mol1->[0]{Name} = $mol2->[0]{Name}\n";
        push @edges, [$mol1->[0]{Name}, $mol2->[0]{Name}];
      }
    }
  }
  clusterization(\@nodes,\@edges);
}
elsif ($renum) {
  write_molden($mols2[0]);
  for (my $i=1; $i<@mols2; $i++) {
    tenzor_sort($mols2[$i-1],$mols2[$i],1e12);
    write_molden($mols2[$i]);
  }
}
elsif ($pair) {
  for (my $i=0; $i<@mols1; $i++) {
    my $mol1 = $mols1[$i];
    my $mol2 = $mols2[$i];
    if (get_formula($mol1) ne get_formula($mol2)) {
      my $ii = $i+1;
      warn "$mol1->[0]{Name} and $mol2->[0]{Name} (pair $ii) are not isomers\n";
      next;
    }
    my $rms;
    if    ($dih)  { $rms = dih_rms($mol1,$mol2) }
    elsif ($dist) { $rms = GrigoryanSpringborg($mol1,$mol2) }
    else          { $rms = tenzor_sort($mol1,$mol2) }
    my $vs = $rms<$RMS ? '==' : '!=';
    if ($dih) {
      printf "$mol1->[0]{Name} $vs $mol2->[0]{Name}: dihedral_diff < %3.0f  %s => %s\n",
             $rms, $mol1->[0]{Energy}||'', $mol2->[0]{Energy}||'';
    }
    elsif ($dist) {
      printf "$mol1->[0]{Name} $vs $mol2->[0]{Name}: diff = %.3f  %s => %s\n",
             $rms, $mol1->[0]{Energy}||'', $mol2->[0]{Energy}||'';
    }
    else {
      printf "$mol1->[0]{Name} $vs $mol2->[0]{Name}: RMSD = %.3f  %s => %s\n",
             $rms, $mol1->[0]{Energy}||'', $mol2->[0]{Energy}||'';
    }
  }
}
else {
  foreach my $mol1 (@mols1) {
    foreach my $mol2 (@mols2) {
      my $rms;
      if    ($dih)  { $rms = dih_rms($mol1,$mol2) }
      elsif ($dist) { $rms = GrigoryanSpringborg($mol1,$mol2) }
      else          { $rms = tenzor_sort($mol1,$mol2) }
      printf "$mol1->[0]{Name} vs. $mol2->[0]{Name}: diff=%.3f\n",$rms if $verbose;
      if ($rms < $RMS) {
        print "$mol1->[0]{Name} = $mol2->[0]{Name}\n";
      }
    }
  } 
}

sub GrigoryanSpringborg {
  my ($mol1,$mol2) = @_;
  my $dist1 = $mol1->[0]{Dist};
  #pp $dist1;
  my $dist2 = $mol2->[0]{Dist};
  #dd $dist1,$dist2;
  return undef if scalar(@$dist1) != scalar(@$dist2);
  my $sum = 0;
  for (my $i=0; $i<@$dist1; $i++) {
    return undef if $dist1->[$i][0] ne $dist2->[$i][0];
    my @di1 = @{$dist1->[$i][1]};
    my @di2 = @{$dist2->[$i][1]};
    my $N = scalar(@di1);
    return undef if $N != scalar(@di2);
    my $sumi = 0;
#    my $aver1 = 1;
#    my $aver2 = 1;
    my $aver1 = $dist1->[$i][3];
    my $aver2 = $dist2->[$i][3];
    for (my $j=0; $j<$N; $j++) {
      $sumi += ($di1[$j]/$aver1 - $di2[$j]/$aver2)**2;
    }
#    if ($N != 1) {
#      my $aver1 = sum(@di1)/$N;
#      my $aver2 = sum(@di2)/$N;
#      for (my $j=0; $j<$N; $j++) {
#        $sumi += ($di1[$j]/$aver1 - $di2[$j]/$aver2)**2;
#      }
#    }
#    else {
#      my ($a1,$a2) = $dist1->[$i][0]=~/(\w+)-(\w+)/;
#      $sumi = (($di1[0]-$di2[0])/($radius{$a1}||2+$radius{$a2}||2)/2)**2;
#    }
    my $m = $dist1->[$i][2];
    $sum += $sumi/$N*$m;
    #printf "%5s (%3d): %6.4f\n", $dist1->[$i][0], $N, sqrt($sumi/$N);
  }
  return sqrt($sum);
}

sub distance_matrix {
 # Take molecules, add sorted distance matrix grouped by bond types
 # Example: trans-diazene H-N=N-H
 # Take array reference
 # [
 #   { Energy => -110.54195368845, Symmetry => "C2h"}, #Properties
 #   ["N", -0.095, -0.621, 0.000, ""],
 #   ["N",  0.095,  0.621, 0.000, ""],
 #   ["H",  0.864, -1.044, 0.000, ""],
 #   ["H", -0.864,  1.044, 0.000, ""],
 # ]
 # Add hash element to Properties
 # Dist => [
 #           ["N-N", [1.256]],
 #           ["N-H", [1.048, 1.048, 1.834, 1.834]],
 #           ["H-H", [2.710]],
 #         ]

  foreach my $mol (@_) {
    #pp $mol;
    my $N = @$mol;
    my %hhh;
    for (my $i=1; $i<$N; $i++) {
      my ($ai,$xi,$yi,$zi) = @{$mol->[$i]};
      for (my $j=$i+1; $j<$N; $j++) {
        my ($aj,$xj,$yj,$zj) = @{$mol->[$j]};
        my $d2 = ($xi-$xj)**2 + ($yi-$yj)**2 + ($zi-$zj)**2;
        #printf "$ai$i-$aj$j\t%6.2f\n", sqrt($d2);
        my $key = $element{$ai}>$element{$aj} ? "$ai-$aj" : "$aj-$ai";
        $key = 'XX-XX' if $W;
        push @{$hhh{$key}}, sqrt($d2);
      }
    }
    
    foreach (keys %hhh) {
      @{$hhh{$_}} = sort {$a <=> $b} @{$hhh{$_}};
      my ($a1,$a2) = m/(\w+)-(\w+)/;
      my $m = 2*$massa{$a1}*$massa{$a2}/($massa{$a1}+$massa{$a2});
      my $aver = sum(@{$hhh{$_}})/@{$hhh{$_}};
      $aver = ($radius{$a1}||2+$radius{$a2}||2)/2 if @{$hhh{$_}}==1;
      $hhh{$_} = [$hhh{$_},$element{$a1},$element{$a2},$m,$aver];
    }
    
    my @dist =
    map {[$_,$hhh{$_}[0],$hhh{$_}[3],$hhh{$_}[4]]} 
    sort {$hhh{$b}[1]<=>$hhh{$a}[1] || $hhh{$b}[2]<=>$hhh{$a}[2]} 
    keys %hhh;
    #pp @dist;
    
    print "=" if $verbose;
    $mol->[0]{Dist} = \@dist;
  }
  print "\n" if $verbose;
}

sub tenzor {
 ###########################################################################
 ## Принимает массив ссылок (молекул) и модифицирует их по месту, 
 ## перемещая в центр масс и ориентируя по осям тензора инерции.
 ## Использует внешний хэш %massa и функцию jacobi().
 ###########################################################################
  my @eigen_num;
  
  foreach my $mol (@_) {
    #dd $mol;
    my (@m,@atom,@x,@y,@z);
    my $N = $#{$mol};
    # Создаем массив атомных масс
    for (my $i=1; $i<=$N; $i++) {
      $atom[$i] = $mol->[$i][0];
      $x[$i] = $mol->[$i][1];
      $y[$i] = $mol->[$i][2];
      $z[$i] = $mol->[$i][3];
      $m[$i] = $massa{$atom[$i]};
    }
    # Вычисляем молекулярную массу
    my $M = sum(@m[1..$N]);
    die "Molecular mass of $mol->[0]{Name} is 0\n" unless $M;
    #print "$M\n";

    # Перемещаем геометрию в центр масс
    my $sum;
    #&write_molden($mol, 'ttt');
    $sum = sum(map { $x[$_]*$m[$_] } 1..$N)/$M;
    $_ -= $sum for @x[1..$N];
    $sum = sum(map { $y[$_]*$m[$_] } 1..$N)/$M;
    $_ -= $sum for @y[1..$N];
    $sum = sum(map { $z[$_]*$m[$_] } 1..$N)/$M;
    $_ -= $sum for @z[1..$N];

    #  Суть: молекулу представляем в виде эллипсоида инерции.
    #  Длина главных осей этого эллипсоида определяются собственными значениями 
    #  тензора инерции, а их направление - соответствующими собственными векторами. 
    #  Г. Голдстейн. Классическая механика. Глава 5.

    # Вычисляем тензор инерции
    # | Ixx Ixy Ixz |
    # | Ixy Iyy Iyz |
    # | Ixz Iyz Izz |
    my (@r2,$Ixx,$Ixy,$Ixz,$Iyy,$Iyz,$Izz);
    $r2[$_] = $x[$_]**2+$y[$_]**2+$z[$_]**2 for 1..$N;
    $Ixx = sum(map {$m[$_]*($r2[$_]-$x[$_]**2)} 1..$N);
    $Iyy = sum(map {$m[$_]*($r2[$_]-$y[$_]**2)} 1..$N);
    $Izz = sum(map {$m[$_]*($r2[$_]-$z[$_]**2)} 1..$N);
    $Ixy = -sum(map {$m[$_]*$x[$_]*$y[$_]} 1..$N);
    $Ixz = -sum(map {$m[$_]*$x[$_]*$z[$_]} 1..$N);
    $Iyz = -sum(map {$m[$_]*$y[$_]*$z[$_]} 1..$N);

    #($Ixx,$Ixy,$Ixz,$Iyy,$Iyz,$Izz) = (6,-2,2,5,0,7);

    #print "Тензор инерции:\n";
    #printf "%12.8f %12.8f %12.8f\n%12.8f %12.8f %12.8f\n%12.8f %12.8f %12.8f\n",
    #         $Ixx,  $Ixy,  $Ixz,   $Ixy,  $Iyy,  $Iyz,   $Ixz,  $Iyz,  $Izz;

    my @eigen = jacobi($Ixx,$Ixy,$Ixz,$Iyy,$Iyz,$Izz) or die "Eigen error\n";
    
    push @eigen_num, [sqrt($eigen[0][0]/$M),
                      sqrt($eigen[1][0]/$M),
                      sqrt($eigen[2][0]/$M)];
    
    # for debug
#    for (my $i=0; $i<=2; $i++) {
#      printf STDERR "%10.6f: %10.4f %10.4f %10.4f\n", $eigen_num[-1][$i], @{$eigen[$i][1]};
#    }
#    print STDERR "\n";
    
    # Преобразуем геометрию к новой системе координат (собственные вектора)
    my @new_geom;
    for (my $i=0; $i<$N; $i++) {
      for (0..2) {
        $new_geom[$_][$i] = $x[$i+1]*$eigen[$_][1][0] + 
                            $y[$i+1]*$eigen[$_][1][1] + 
                            $z[$i+1]*$eigen[$_][1][2];
      }
    }
    #write_molden();

    # Переставляем оси эллипсоида, чтобы хорошо смотрелось в молдене
    @x[1..$N] = @{$new_geom[1]};
    @y[1..$N] = @{$new_geom[2]};
    @z[1..$N] = @{$new_geom[0]};
  
    # Из новой геометрии делаем две молекулы (одну с отражением по z)
    my ($mol1,$mol2);
    for (my $i=1; $i<=$N; $i++) {
      $mol2->[$i][0] = $mol1->[$i][0] = $mol->[$i][0];
      $mol2->[$i][1] = $mol1->[$i][1] = $x[$i];
      $mol2->[$i][2] = $mol1->[$i][2] = $y[$i];
      $mol1->[$i][3] =  $z[$i];
      $mol2->[$i][3] = -$z[$i];
    }
    
    my $znakz = 1;
    #my $rms1 = get_rms($mol,$mol1);
    my $rms2 = get_rms($mol,$mol2);
    print '=' if $verbose;
    #printf "%20.10f %20.10f\n", $rms1, $rms2;
    $znakz = -1 if $rms2 < 1e-8;
  
    # Перезаписываем молекулу
    for (my $i=1; $i<=$N; $i++) {
      $mol->[$i][1] = $x[$i];
      $mol->[$i][2] = $y[$i];
      $mol->[$i][3] = $znakz*$z[$i];
    }
    #dd $mol;
    
    # Проверка на вырожденность собственных чисел тензора инерции
    # ... по плоскости xy
    #dd @eigen_num;
    if ($eigen_num[-1][1]-$eigen_num[-1][2] < 0.001) {
      #dd $mol;
      my @d = grep {$_->[0] > 0.1} proections($mol,'xy');
      #dd @d;
      # Если есть хотя бы 3 равноудаленных атома...
      if (@d>2 && $d[1][0]-$d[0][0]<0.001 && $d[2][0]-$d[0][0]<0.001) {
        # будем переориентировать, направив ось x к первому из них
        my $n = $d[0][1];
        my $cos = $mol->[$n][1]/sqrt($mol->[$n][1]**2+$mol->[$n][2]**2);
        my $sin = sqrt(1-$cos**2);
        $sin = -$sin if $mol->[$n][2] < 0;
        for (my $i=1; $i<@$mol; $i++) {
          my $x =  $cos*$mol->[$i][1] + $sin*$mol->[$i][2];
          my $y = -$sin*$mol->[$i][1] + $cos*$mol->[$i][2];
          $mol->[$i][1] = $x;
          $mol->[$i][2] = $y;
        }
      }
    }
    # ... по плоскости xz
    if ($eigen_num[-1][0]-$eigen_num[-1][1] < 0.001) {
      #dd $mol;
      my @d = grep {$_->[0] > 0.1} proections($mol,'xz');
      #dd @d;
      # Если есть хотя бы 3 равноудаленных атома...
      if (@d>2 && $d[1][0]-$d[0][0]<0.001 && $d[2][0]-$d[0][0]<0.001) {
        # будем переориентировать, направив ось x к первому из них
        my $n = $d[0][1];
        my $cos = $mol->[$n][1]/sqrt($mol->[$n][1]**2+$mol->[$n][3]**2);
        my $sin = sqrt(1-$cos**2);
        $sin = -$sin if $mol->[$n][2] < 0;
        for (my $i=1; $i<@$mol; $i++) {
          my $x =  $cos*$mol->[$i][1] + $sin*$mol->[$i][3];
          my $z = -$sin*$mol->[$i][1] + $cos*$mol->[$i][3];
          $mol->[$i][1] = $x;
          $mol->[$i][3] = $z;
        }
      }
      #dd $mol;
    }
  }
  print "\n" if $verbose;
  return @eigen_num;
}
  
sub tenzor_sort {
 ###########################################################################
 ## Принимает две ссылки (молекулы) и sort_RMS третьим параметром.
 ## Молекулы должны быть предварительно пропущены через tenzor.
 ## Возвращает rms между этими молекулами.
 ## Если rms < sort_RMS, то порядок атомов второй молекулы сортируется по первой.
 ## rms = sqrt(sum(m*r**2)/sum(m))
 ## Использует функции znak() и kinen()
 ###########################################################################
  my ($mol0, $mol, $sort_RMS) = @_;
  $sort_RMS ||= 0;
  my $N = $#{$mol0};
  # Создаем структуру данных "молекула", 
  # объединяющую символы атомов, их порядок в геометрии и координаты.
  # Это хэш (атомных символов) хэшей (их номера от 1 до N) массивов (координат)
  my (%mol0, %mol);
  foreach (1..$N) {
    my @a = @{$mol0->[$_]};
    $mol0{$a[0]}->{$_} = [@a[1..$#a]];
    @a = @{$mol->[$_]};
    $mol{$a[0]}->{$_} = [@a[1..$#a]];
  }
  
  # Нужно согласовать направления осей текущей
  # геометрии с направлениями предыдущей. Эти направления после уже 
  # проделанных преобразования могут отличаться только знаками. Сложность
  # же заключается в том, что порядок нумерации атомов в двух геометриях
  # может быть различен. Но предполагаем, что две соседние геометрии 
  # не сильно отличаются друг от друга
  
  # Вычисляем "кинетич. энергию" между текущей и предыдущей молекулами
  # В возвращаемом массиве 0-й элемент - это "кинетич. энергия", остальные - 
  # пермутационный вектор, определяющий порядок нумерации текущей молекулы,
  # при котором она лучше всего накладывается на предыдущую.
  my @jmin = kinen(\%mol,\%mol0);

  # Далее будем последовательно менять знаки координат текущей молекулы
  # (отражение) и смотреть, при каком отражении будет мин. кин. энергия
  # Последовательность отражений:
  #  0 +x +y +z    
  #  1 +x -y +z    
  #  2 -x -y +z    
  #  3 -x +y +z    
  #  4 -x +y -z    
  #  5 +x +y -z    
  #  6 +x -y -z    
  #  7 -x -y -z    

  my $reflect = 0;    # Счетчик отражений
  my $reflectmin = 0; # Номер отражения с минимальн. кинетич. энергией

  # Для последовательности отражений, подобранной так, чтобы меняя каждый раз 
  # знак только одной координаты (0 - x, 1 - y, 2 - z), охватить все варианты
  my @j;
  foreach (1,0,1,2,0,1,0) {
    # Меняем знак соответствующей координаты
    znak(\%mol, $_);
    $reflect++;  # Счетчик отражений
    # Если нужно сохранить конфигурацию, пропускаем вычисления на нечетных шагах
    # (останутся только повороты вокруг осей z,x,y)
    #  0 +x +y +z    
    #  2 -x -y +z    
    #  4 -x +y -z    
    #  6 +x -y -z    
    if ($chiral) {
      next if $reflect%2;
    }
    # Вычисляем кин. энергию
    @j = kinen(\%mol,\%mol0);
    #print "@j\n";
    # Смотрим, не лучше ли очередная кин. энергия
    if ($j[0] < $jmin[0]) {
      @jmin = @j;
      $reflectmin = $reflect;
    }
  }
  #warn "$reflectmin @jmin\n";

  # Вспомогательный массив, определяющий, как нужно поменять знаки координат,
  # чтобы из молекулы, полученной после последнего отражения (-x -y -z), 
  # восстановить молекулу после i-го отражения
  my @ar = ([0,1,2],[0,2],[2],[1,2],[1],[0,1],[0],[]);

  # Меняем знаки, чтобы получить наименьшую кин. энергию 
  znak(\%mol, @{$ar[$reflectmin]});

  # Переносим полученные координаты из "молекулы" %mol в массивы @x,@y,@z
  # в зависимости от значения опции -sort
  if ($sort_RMS > $jmin[0]) {       # Сортировка по предыдущей геометрии
    my $i = 1;
    #print "@jmin\n";
    foreach my $j (@jmin[1..$#jmin]) {
      foreach my $atom (keys %mol) {
        if (exists $mol{$atom}{$j}) {
          @{$mol->[$i]} = ($atom, @{$mol{$atom}{$j}});
          $i++;
        }
      }
    }
  }
  print "Permutation vector:", join(',', @jmin[1..$#jmin]), "\n" if $verbose;
  #write_molden($mol0,$mol);
  return $jmin[0];  
}

sub znak {
 ###########################################################################
 ## Менять знаки координат в структуре $h{$atom}{$index}[$xyz]
 ## Usage: znak(\%h, 0|1|2); x - 0, y - 1, z - 2
 ###########################################################################
  my %h = %{$_[0]} and shift;
  foreach my $xyz (@_) {
    foreach my $atom (keys %h) {
      for (keys %{$h{$atom}}) {
        $h{$atom}{$_}[$xyz] *= -1;
      }
    }
  }
}

sub kinen {
 ###########################################################################
 ## Вычисляет "кинетич. энергию" между сравниваемой и предыдущей молекулами
 ## Принимает две ссылки на сравниваемую и предыдущую молекулы. Молекула -
 ## это хэш (атомных символов) хэшей (их номера от 1 до N) массивов (координат)
 ## В возвращаемом массиве 0-й элемент - это "кинетич. энергия", остальные - 
 ## пермутационный вектор, определяющий порядок нумерации текущей молекулы,
 ## при котором она лучше всего накладывается на предыдущую.
 ## "кинетич. энергия" rms = sqrt(sum(m*r**2)/sum(m))
 ## Использует внешний хэш %massa
 ###########################################################################

  # Делаем независимую копию 1-й (сравниваемой) молекулы
  my %mol;
  foreach my $atom (keys %{$_[0]}) {
    while (my ($ind,$aref) = each %{$_[0]{$atom}}) {
      $mol{$atom}{$ind} = [@{$aref}];
    }
  }
  
  # 2-я (предыдущая) молекула
  my %mol0 = %{$_[1]};
  
  # Вычисляем молекулярную массу
  my $M = sum(map {$massa{$_}*(keys %{$mol{$_}})} keys %mol);
  
  my $E = 0;
  my @j;
  
  # Для углеродов, водородов и т.д. начиная с более тяжелых
  foreach my $atom (sort {$massa{$b}<=>$massa{$a}} keys %mol0) {
    
    # Если в молекулах разное кол-во, например, углеродов, то выход
    if (scalar(keys %{$mol{$atom}}) != scalar(keys %{$mol0{$atom}})) {
      warn "Not isomers ($atom)\n";
      return;
    }
    
    # Для каждого из атомов данного сорта 2-й молекулы
    # отсортированных по увеличению расстояния до центра масс
    foreach my $i (map {$_->[0]}
                   sort {$a->[1]<=>$b->[1]} 
                   map {[$_,$mol0{$atom}{$_}[0]**2+$mol0{$atom}{$_}[1]**2+$mol0{$atom}{$_}[2]**2]}
                   keys %{$mol0{$atom}}) {
      my $rmin = 1e10;
      my $jmin;
      # перебираем атомы того же сорта 1-й молекулы
      foreach my $j (keys %{$mol{$atom}}) {
        # вычисляем расстояние между атомами двух молекул
        #my $r = sum(map {($mol{$atom}{$j}[$_] - $mol0{$atom}{$i}[$_])**2} 0..2);
        my $r = ($mol{$atom}{$j}[0] - $mol0{$atom}{$i}[0])**2 +
                ($mol{$atom}{$j}[1] - $mol0{$atom}{$i}[1])**2 +
                ($mol{$atom}{$j}[2] - $mol0{$atom}{$i}[2])**2;
          
        # запоминая ближайший атом 1-й молекулы
        if ($r<$rmin) {
          $rmin = $r;
          $jmin = $j;
        }
      }
      # Сумимируем кинетическую энергию
      $E += $rmin*$massa{$atom};
      # i-му атому 2-й молекулы ставим в соответствие ближайший атом из 1-й
      $j[$i] = $jmin;
      # и удаляем его из 1-й молекулы
      delete $mol{$atom}{$jmin};
    }
    # Записываем накопленную кин. энергию в начало массива
    $j[0] = $E;
  }
  $j[0] = sqrt($j[0]/$M);
  # Возвращаем массив
  return @j;
}

sub jacobi {
 ###########################################################################
 ## Собственные числа и векторы симметричной матрицы методом Якоби.
 ## Код из Жориной программы dte (C).
 ## Принимает массив ($Ixx,$Ixy,$Ixz,$Iyy,$Iyz,$Izz).
 ## Возвращает отсортированный (по убыванию собств. чисел) массив двухэлементных
 ## массивов. 1-й элемент - собств. число, 2-й - соответствующий ему вектор.
 ###########################################################################
  my @a = ([ $_[0],$_[1],$_[2] ],
           [ $_[1],$_[3],$_[4] ],
           [ $_[2],$_[4],$_[5] ]);
  my ($j, $iq, $ip, $i, $n);
  my ($tresh, $theta, $tau, $t, $sm, $s, $h, $g, $c);
  my $maxit = 50;

  $n = 3;

  my @v = ([1,0,0], [0,1,0], [0,0,1]);
  my @b = my @d = ($a[0][0],$a[1][1],$a[2][2]);
  my @z = (0,0,0);
  
  my $nrot = 0;
  for ($i=0; $i<$maxit; $i++) {
    
    $sm = abs($a[0][1]) + abs($a[0][2]) + abs($a[1][2]);
    
    if ($sm == 0) {
      return sort {$b->[0] <=> $a->[0]}
              ([$d[0], [$v[0][0],$v[1][0],$v[2][0]]],
              [$d[1], [$v[0][1],$v[1][1],$v[2][1]]],
              [$d[2], [$v[0][2],$v[1][2],$v[2][2]]]);
    }
    # The normal return, which relies on 
    # quadratic convergence to machine underflow. 
    
    $tresh = ($i < 4) ? 0.2*$sm/($n*$n) : 0;
    # ...on the first three sweeps. 
    # ...thereafter. 
    
    for ($ip=0; $ip<$n-1; $ip++) {
      for ($iq=$ip+1; $iq<$n; $iq++) {
        $g = 100*abs($a[$ip][$iq]);
        # After four sweeps, skip the rotation if the off-diagonal element is small. 
        if ($i>4 && abs($d[$ip])+$g == abs($d[$ip]) && abs($d[$iq])+$g == abs($d[$iq])) {
          $a[$ip][$iq] = 0;
        }
        elsif (abs ($a[$ip][$iq]) > $tresh) {
          $h = $d[$iq]-$d[$ip];
          if (abs($h)+$g == abs($h)) {
            $t = $a[$ip][$iq]/$h;
          }  # $t = 1/(2theta) 
          else
          {
            $theta = 0.5*$h/$a[$ip][$iq];
            $t = 1/(abs($theta)+sqrt(1+$theta*$theta));
            if ($theta < 0) {$t = -$t;}
          }
          $c = 1/sqrt(1+$t*$t);
          $s = $t*$c;
          $tau = $s/(1+$c);
          $h = $t*$a[$ip][$iq];
          $z[$ip] -= $h;
          $z[$iq] += $h;
          $d[$ip] -= $h;
          $d[$iq] += $h;
          $a[$ip][$iq] = 0;
          for ($j=0; $j<$ip; $j++) {  # Case of rotations 0 <= $j < p. 
            $g = $a[$j][$ip]; 
            $h = $a[$j][$iq];
            $a[$j][$ip] = $g-$s*($h+$g*$tau);
            $a[$j][$iq] = $h+$s*($g-$h*$tau);
          }
          for ($j=$ip+1; $j<$iq; $j++) {  # Case of rotations p < $j < q. 
            $g = $a[$ip][$j]; 
            $h = $a[$j][$iq];
            $a[$ip][$j] = $g-$s*($h+$g*$tau);
            $a[$j][$iq] = $h+$s*($g-$h*$tau);
          }
          for ($j=$iq+1; $j<$n; $j++) {  # Case of rotations q < $j < $n. 
            $g = $a[$ip][$j]; 
            $h = $a[$iq][$j];
            $a[$ip][$j] = $g-$s*($h+$g*$tau);
            $a[$iq][$j] = $h+$s*($g-$h*$tau);
          }
          for ($j=0; $j<$n; $j++) {
            $g = $v[$j][$ip]; 
            $h = $v[$j][$iq];
            $v[$j][$ip] = $g-$s*($h+$g*$tau);
            $v[$j][$iq] = $h+$s*($g-$h*$tau);
          }
          ++$nrot;
        }
      }
    }
    for ($ip=0; $ip<$n; $ip++) {
      $b[$ip] += $z[$ip];
      $d[$ip] = $b[$ip];      # Update $d with the sum of tapq, 
      $z[$ip] = 0;        # and reinitialize $z. 
    }
  }
  return undef;
}

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;
}

sub write_molden {
  my $oldfh;
  my ($full_title,$ppm);
  
  if (@_ > 1 && !ref($_[-1])) {
    my $file = pop @_;
    open F, '>', $file or do {warn "Can't write to $file: $!\n"; return};
    $oldfh = select F;
  }
  
  foreach my $mol (@_) {
    #pp $mol;
    my $N = $#{$mol};
    print " $N\n";
	  if ($full_title) {
      print "$mol->[0]{Title}";
    }
    else {
      foreach my $f (qw/Energy Charge Mult Symmetry Dipole HoF Edisp Ellips sigma Content ZPE g-tenzor/) {
        #print "  $f $mol->[0]{$f}";
        print "  $f $mol->[0]{$f}" if defined $mol->[0]{$f};
      } 
      if (defined $mol->[0]{G}) {
        foreach my $g (@{$mol->[0]{G}}) {
          next unless $g->[0];
          print "  G($g->[0]) $g->[1]";
        } 
      }
      print "\n";
    }
	  for (my $i=1; $i<=$N; $i++) {
      printf " %-2s %12.8f %12.8f %12.8f", @{$mol->[$i]}[0..3];
      print "    $mol->[$i][4]" if $ppm && $mol->[$i][4];
      print "\n";
	  }
  }
  
  if ($oldfh) {
    close F;
    select $oldfh;
  }
}

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

sub get_rms {
 ############################################################################
 ## Переделанная на перл coord.c Г.Сальникова
 ## Принимает две ссылки на геометрии (см. read_molden() из conformers)
 ## Возвращает RMSD между геометриями
 ## Во вторую ссылку записывает новую, сглаженную геометрию
 ## Использует глобальный хэш %massa с массами элементов (см. smooth)
 ############################################################################
  my ($mol0, $mol1) = @_;
  my ($i, $n, $rot);
  my (@m, @x0, @y0, @z0, @x1, @y1, @z1, $xc, $yc, $zc, $e, $e0,
    $mtot, $tg1, $tg2, $phi, $phix, $phiy, $phiz);
  my $M_PI = atan2(0,-1);
  my $FLT_EPSILON = 1e-8;
  my $DEBUG = 0;
  
  #/*** Input data preparation ***/
  $n = $#{$mol0};
  for ($i=0; $i<$n; $i++) {
    $m[$i] = $massa{$mol0->[$i+1][0]};
    die "Inconsistent data\n" if $massa{$mol1->[$i+1][0]} != $m[$i];
    $mtot += $m[$i];
    $x0[$i] = $mol0->[$i+1][1];
    $y0[$i] = $mol0->[$i+1][2];
    $z0[$i] = $mol0->[$i+1][3];
    $x1[$i] = $mol1->[$i+1][1];
    $y1[$i] = $mol1->[$i+1][2];
    $z1[$i] = $mol1->[$i+1][3];
  }
  
  #/*** 1-st translation to center of mass ***/

  $xc = $yc = $zc = 0;
  for ($i=0; $i<$n; $i++)
  {
    $xc += $m[$i]*$x0[$i];
    $yc += $m[$i]*$y0[$i];
    $zc += $m[$i]*$z0[$i];
  }
  $xc /= $mtot;
  $yc /= $mtot;
  $zc /= $mtot;

  for ($i=0; $i<$n; $i++)
  {
    $x0[$i] -= $xc;
    $y0[$i] -= $yc;
    $z0[$i] -= $zc;
  }

  if ($DEBUG) {
    printf ("1-st molecule in center of mass\n");
    for ($i=0; $i<$n; $i++) {
      printf "%2.0f   %10f   %10f   %10f\n", $m[$i], 
              $x0[$i], $y0[$i], $z0[$i];
    }
    printf "1-st center of mass translation on: (%10f, %10f, %10f)\n",
            -$xc, -$yc, -$zc;
  }

  #/*** 2-nd translation to center of mass ***/

  $xc = $yc = $zc = 0;
  for ($i=0; $i<$n; $i++)
  {
    $xc += $m[$i]*$x1[$i];
    $yc += $m[$i]*$y1[$i];
    $zc += $m[$i]*$z1[$i];
  }
  $xc /= $mtot;
  $yc /= $mtot;
  $zc /= $mtot;

  for ($i=0; $i<$n; $i++)
  {
    $x1[$i] -= $xc;
    $y1[$i] -= $yc;
    $z1[$i] -= $zc;
  }

  if ($DEBUG) {
    printf ("2-nd molecule in center of mass\n");
    for ($i=0; $i<$n; $i++) {
      printf "%2.0f   %10f   %10f   %10f\n", 
              $m[$i], $x1[$i], $y1[$i], $z1[$i];
    }
    printf "2-nd center of mass translation on: (%10f, %10f, %10f)\n",
            -$xc, -$yc, -$zc;
  }

  $e = 0;
  for ($i=0; $i<$n; $i++) {
    $e += $m[$i]*(($x1[$i]-$x0[$i])*($x1[$i]-$x0[$i])+
         ($y1[$i]-$y0[$i])*($y1[$i]-$y0[$i])+
         ($z1[$i]-$z0[$i])*($z1[$i]-$z0[$i]));
  }
  printf ("2*energy: %g\n", $e) if $DEBUG;

  for ($rot=1; ; $rot++)
  {
    $e0 = $e;

    #/*** Rotation around X ***/

    $tg1 = $tg2 = 0;
    for ($i=0; $i<$n; $i++)
    {
      $tg1 += $m[$i]*($y0[$i]*$z1[$i]-$z0[$i]*$y1[$i]);
      $tg2 += $m[$i]*($y0[$i]*$y1[$i]+$z0[$i]*$z1[$i]);
    }
    $phi = atan2 ($tg1, $tg2);

    for ($i=0; $i<$n; $i++)
    {
      $yc = $y1[$i]*cos($phi)+$z1[$i]*sin($phi);
      $zc = $z1[$i]*cos($phi)-$y1[$i]*sin($phi);
      $y1[$i] = $yc;
      $z1[$i] = $zc;
    }

    $e = 0;
    for ($i=0; $i<$n; $i++) {
      $e += $m[$i]*(($x1[$i]-$x0[$i])*($x1[$i]-$x0[$i])+
     ($y1[$i]-$y0[$i])*($y1[$i]-$y0[$i])+
     ($z1[$i]-$z0[$i])*($z1[$i]-$z0[$i]));
    }
    if ($DEBUG) {
      for ($i=0; $i<$n; $i++) {
        printf "%2.0f   %10f   %10f   %10f\n", 
                $m[$i], $x1[$i], $y1[$i], $z1[$i];
      }
      printf "after %d-th rotation around X on: %g (%g deg)\n2*energy: %g\n",
              $rot, $phi, $phi*180/$M_PI, $e;
    }

    $phix = $phi;

    #/*** Rotation around Y ***/

    $tg1 = $tg2 = 0;
    for ($i=0; $i<$n; $i++)
    {
      $tg1 += $m[$i]*($z0[$i]*$x1[$i]-$x0[$i]*$z1[$i]);
      $tg2 += $m[$i]*($z0[$i]*$z1[$i]+$x0[$i]*$x1[$i]);
    }
    $phi = atan2 ($tg1, $tg2);

    for ($i=0; $i<$n; $i++)
    {
      $xc = $x1[$i]*cos($phi)-$z1[$i]*sin($phi);
      $zc = $z1[$i]*cos($phi)+$x1[$i]*sin($phi);
      $x1[$i] = $xc;
      $z1[$i] = $zc;
    }

    $e = 0;
    for ($i=0; $i<$n; $i++) {
      $e += $m[$i]*(($x1[$i]-$x0[$i])*($x1[$i]-$x0[$i])+
     ($y1[$i]-$y0[$i])*($y1[$i]-$y0[$i])+
     ($z1[$i]-$z0[$i])*($z1[$i]-$z0[$i]));
    }
    if ($DEBUG) {
      for ($i=0; $i<$n; $i++) {
        printf "%2.0f   %10f   %10f   %10f\n", 
                $m[$i], $x1[$i], $y1[$i], $z1[$i];
      }
      printf "after %d-th rotation around Y on: %g (%g deg)\n2*energy: %g\n",
              $rot, $phi, $phi*180/$M_PI, $e;
    }

    $phiy = $phi;

    #/*** Rotation around Z ***/

    $tg1 = $tg2 = 0;
    for ($i=0; $i<$n; $i++)
    {
      $tg1 += $m[$i]*($x0[$i]*$y1[$i]-$y0[$i]*$x1[$i]);
      $tg2 += $m[$i]*($x0[$i]*$x1[$i]+$y0[$i]*$y1[$i]);
    }
    $phi = atan2 ($tg1, $tg2);

    for ($i=0; $i<$n; $i++)
    {
      $xc = $x1[$i]*cos($phi)+$y1[$i]*sin($phi);
      $yc = $y1[$i]*cos($phi)-$x1[$i]*sin($phi);
      $x1[$i] = $xc;
      $y1[$i] = $yc;
    }

    $e = 0;
    for ($i=0; $i<$n; $i++) {
      $e += $m[$i]*(($x1[$i]-$x0[$i])*($x1[$i]-$x0[$i])+
     ($y1[$i]-$y0[$i])*($y1[$i]-$y0[$i])+
     ($z1[$i]-$z0[$i])*($z1[$i]-$z0[$i]));
    }
    if ($DEBUG) {
      for ($i=0; $i<$n; $i++) {
        printf "%2.0f   %10f   %10f   %10f\n", 
               $m[$i], $x1[$i], $y1[$i], $z1[$i];
      }
      printf "after %d-th rotation around Z on: %g (%g deg)\n2*energy: %g\n",
        $rot, $phi, $phi*180/$M_PI, $e;
    }

    $phiz = $phi;

    last if $rot > 100;
    last if (abs($phix) < $FLT_EPSILON &&
             abs($phiy) < $FLT_EPSILON &&
             abs($phiz) < $FLT_EPSILON &&
             ($e == $e0 || 
             abs($e+$e0) < $FLT_EPSILON || 
             (abs(($e-$e0)/($e+$e0)) < $FLT_EPSILON)));
    if ($DEBUG) {
      print "
last if (abs($phix) < $FLT_EPSILON &&
         abs($phiy) < $FLT_EPSILON &&
         abs($phiz) < $FLT_EPSILON &&
         ($e == $e0 || 
         abs($e+$e0) < $FLT_EPSILON || 
         (abs(($e-$e0)/($e+$e0)) < $FLT_EPSILON)));
      \n";
    }
  }

  #/*** Output data preparation ***/

  for ($i=0; $i<$n; $i++) {
    $mol1->[$i+1][1] = $x1[$i];
    $mol1->[$i+1][2] = $y1[$i];
    $mol1->[$i+1][3] = $z1[$i];
  }
  if ($DEBUG) {
    for ($i=0; $i<$n; $i++) {
      printf "%2.0f   %10f   %10f   %10f\n", $m[$i], $x1[$i], $y1[$i], $z1[$i];
    }
    printf "2*energy: %g\n", $e;
  }

  #/*** Finished ***/
  
#  return sqrt($e);
  return $e;
}

sub delete_eq_end {
  my $aref = shift;
  return '' if @$aref == 1;
  my $end = '';
  while (1) {
    my @a = @$aref;
    my $a0 = substr $a[0], -1, 1, '';
    for (my $i=1; $i<@a; $i++) {
      my $ai = substr $a[$i], -1, 1, '';
      return reverse($end) if $ai ne $a0;
    }
    $end .= $a0;
    @$aref = @a;
  }
}

sub proections {
 ############################################################################
 ## Принимает ссылкy на молекулу и вторым параметром - плоскость, 
 ## на которую делать проекцию ('xy','xz' или 'yz').
 ## Возвращает сортированный по длинам проекций массив 
 ## [проекция на плоскость xy вектора 0,0,0-атом , номер атома]
 ############################################################################
  my $mol = shift;
  my $plane = shift;
  $plane ||= 'xy';
  my @d;
  for (my $i=1; $i<@$mol; $i++) {
    my $d;
    if ($plane eq 'xy') {
      $d = sqrt($mol->[$i][1]**2 + $mol->[$i][2]**2);
    }
    elsif ($plane eq 'xz') {
      $d = sqrt($mol->[$i][1]**2 + $mol->[$i][3]**2);
    }
    elsif ($plane eq 'yz') {
      $d = sqrt($mol->[$i][2]**2 + $mol->[$i][3]**2);
    }
    push @d, [$d,$i];
  }
  @d = sort {$a->[0] <=> $b->[0]} @d;
  return @d;
}

sub get_formula {
	my @mol = @{$_[0]};
	my %f;
	my $formula = '';
	$f{ucfirst lc $_->[0]}++ foreach @mol[1..$#mol];
	foreach (sort by_Hill keys %f) {
	  $formula .= $_ . ($f{$_}==1 ? '' : $f{$_});
	}
	return $formula;
}

sub by_Hill {
	return -1 if uc($a) eq 'C' and uc($b) ne 'C';
	return  1 if uc($a) ne 'C' and uc($b) eq 'C';
	return -1 if uc($a) eq 'H' and uc($b) ne 'H';
	return  1 if uc($a) ne 'H' and uc($b) eq 'H';
	uc($a) cmp uc($b);
}

sub dih_rms {
 # $mol->[0]{Dih} = [d, d1, ...]
  my ($mol1,$mol2) = @_;
  #my $rms2;
  my $max = 0;
  my @dih1 = @{$mol1->[0]{Dih}};
  my @dih2 = @{$mol2->[0]{Dih}};
  die "Inconsistent dihedral lists\n" if @dih1 != @dih2;
  for (my $i=0; $i<@dih1; $i++) {
    my $diff =  $dih1[$i] - $dih2[$i];
    $diff += 360 if $diff <= -180;
    $diff -= 360 if $diff > 180;
    $max = $diff if abs($max)<abs($diff);
    if ($verbose) {
      printf "Dihedral %d: %.0f - %.0f = %.0f\n", $i+1,$dih1[$i],$dih2[$i],$diff;
    }
    #$rms2 += $diff**2;
  }
  close L;
  #return sqrt($rms2/@dih1);
  return abs($max);
}

sub dihedral {
  my ($mol,$d) = @_;
  #d |i-j-k-l| диэдральный угол i-j-k-l (-180..180)|
  my $radian = 57.2957795130823;
  my ($i,$j,$k,$l) = split /,/, $d;
  my ($xi,$yi,$zi) = @{$mol->[$i]}[1..3];
  my ($xj,$yj,$zj) = @{$mol->[$j]}[1..3];
  my ($xk,$yk,$zk) = @{$mol->[$k]}[1..3];
  my ($xl,$yl,$zl) = @{$mol->[$l]}[1..3];
  my $A1 = ($yj-$yi)*($zk-$zi)-($zj-$zi)*($yk-$yi);
  my $B1 = ($zj-$zi)*($xk-$xi)-($xj-$xi)*($zk-$zi);
  my $C1 = ($xj-$xi)*($yk-$yi)-($yj-$yi)*($xk-$xi);
  my $A2 = ($yk-$yj)*($zl-$zj)-($zk-$zj)*($yl-$yj);
  my $B2 = ($zk-$zj)*($xl-$xj)-($xk-$xj)*($zl-$zj);
  my $C2 = ($xk-$xj)*($yl-$yj)-($yk-$yj)*($xl-$xj);
  my $Vec = ($B1*$C2-$C1*$B2)*($xk-$xj)+
         ($C1*$A2-$A1*$C2)*($yk-$yj)+
         ($A1*$B2-$B1*$A2)*($zk-$zj);
  my $sign = ($Vec>=0) ? 1 : -1;
  my $cos = ($A1*$A2+$B1*$B2+$C1*$C2)/sqrt(($A1**2+$B1**2+$C1**2)*($A2**2+$B2**2+$C2**2));
  $cos = 1 if $cos>1 && $cos<1.0000001;
  #return $sign*rad2deg(acos($cos));
  my $acos = atan2(sqrt(1-$cos**2),$cos);
  return $sign*$radian*$acos;
}

sub dihedral2 {
  my ($mol,$d) = @_;
  my $radian = 57.2957795130823;
  my ($i,$j,$k,$l) = split /,/, $d;
  my ($x1,$y1,$z1) = @{$mol->[$i]}[1..3];
  my ($x2,$y2,$z2) = @{$mol->[$j]}[1..3];
  my ($x3,$y3,$z3) = @{$mol->[$k]}[1..3];
  my ($x4,$y4,$z4) = @{$mol->[$l]}[1..3];

  # b1
  my @b_a;
  $b_a[0] = -($x1 - $x2);
  $b_a[1] = -($y1 - $y2);
  $b_a[2] = -($z1 - $z2);

  # b2
  my @b_c;
  $b_c[0] = $x2 - $x3;
  $b_c[1] = $y2 - $y3;
  $b_c[2] = $z2 - $z3;

  # b3
  my @c_d;
  $c_d[0] = $x4 - $x3;
  $c_d[1] = $y4 - $y3;
  $c_d[2] = $z4 - $z3;

  VectorNormalisation(@b_c);
  VectorNormalisation(@b_a);
  VectorNormalisation(@c_d);

  my @n1 = CrossProduct(@b_a, @b_c);
  my @n2 = CrossProduct(@b_c, @c_d);
  my @m = CrossProduct(@n1, @b_c);

  my $x = DotProduct(@n1, @n2);
  my $y = DotProduct(@m, @n2);

  return -$radian*atan2($y, $x);
}

sub VectorNormalisation { 
  my $lenght = sqrt($_[0]**2 + $_[1]**2 + $_[2]**2); 
  $_[0] /= $lenght; 
  $_[1] /= $lenght; 
  $_[2] /= $lenght;
}

sub DotProduct {
  my @v = @_[0..2];
  my @w = @_[3..5];
  return $v[0] * $w[0] + $v[1] * $w[1] + $v[2] * $w[2];
}

sub CrossProduct {
  my @v = @_[0..2];
  my @w = @_[3..5];
  my @cross;
  $cross[0] = $w[1] * $v[2] - $w[2] * $v[1];
  $cross[1] = $w[2] * $v[0] - $w[0] * $v[2];
  $cross[2] = $w[0] * $v[1] - $w[1] * $v[0];
  return @cross;
}


sub calculate_dihedrals {
 # Add $mol->[0]{Dih} = [d, d1, ...]
  my $dih = pop @_;
  my @dih;
  foreach (split ' ', $dih) {
    while (s/(\d+,(\d+,\d+,\d+))/$2/) {
      push @dih, $1;
    }
  }
  foreach my $mol (@_) {
    foreach (@dih) {
      push @{$mol->[0]{Dih}}, dihedral($mol,$_);
    }
  }
}

sub clusterization {
  my @nodes = @{$_[0]};
  my @edges = @{$_[1]};
  my @clusters;
  
  foreach my $edge (@edges) {
    #dd $edge;
    my $found;
    foreach (@clusters) {
      if (exists($_->{$edge->[0]}) || exists($_->{$edge->[1]})) {
        $_->{$edge->[0]}++;
        $_->{$edge->[1]}++;
        $found = 1;
      }
    }
    if (! $found) {
      push @clusters, {$edge->[0]=>1,$edge->[1]=>1};
    }
    #dd @clusters;
  }
  
  # Add non-edged nodes
  my %nodes_in_clusters = map {%$_} @clusters;
  foreach (@nodes) {
    next if exists $nodes_in_clusters{$_};
    push @clusters, {$_=>0};
  }
  
  # Sort @clusters by @ARGV (for output)
  my %nodes = map {$nodes[$_-1]=>$_} 1..@nodes;
  #dd %nodes; exit;
  my @sort_clusters;
  foreach (@clusters) {
    push @sort_clusters, [sort {$nodes{$a}<=>$nodes{$b}} keys %$_];
  }
  @sort_clusters = sort {$nodes{$a->[0]}<=>$nodes{$b->[0]}} @sort_clusters;
  # Output
  print "\nClusterization\n";
  foreach (@sort_clusters) {
    print join(',', @$_), "\n";
  }
}