#!/usr/bin/perl -ws

use strict;
#use Data::Dump 'pp';
our ($h,$help,$ijk,$IJK,$rm_in_core,$rm_dist,$R,$ii,$rot,$make_R); # Опции

if ($h or $help) {
  (my $prog = $0) =~ s/.+\///;
  print <<HELP;
Совмещает молекулы и выбрасывает совпадающие атомы.
Usage: $prog [-ijk=i,j,k] [-IJK=I,J,K] core.xyz substitutent.xyz [> core_subst.xyz]
Зависимости: Perl

Совмещение происходит по тройкам атомов из опций -ijk и -IJK.
По умолчанию -ijk=1,2,3 и -IJK=1,2,3.
-ijk относится к замещаемой молекуле/молекулам (core.xyz), 
-IJK - к заместителю (substitutent.xyz или последняя из структур). 

При совмещении происходит переориентация:
      x     y     z
 i:  0.00  0.00  0.00    Атом i в центре координат
 j:        0.00  0.00    Атом j на оси x
 k:              0.00    Атом k в плоскости xy

Опции:
-ijk=i,j,k  номера тройки атомов в замещаемых молекулах
-IJK=I,J,K  номера тройки атомов в заместителе (последней структуре)
-rm_in_core  накладывающиеся атомы удаляются в остове, а не в заместителе
-rm_dist=0.8 расстояние-признак накладывающихся атомов (ангстремы)

Использование в режиме аннелирования:
$prog -ijk=2,3,12 -IJK=1,2,3 Indole.xyz Benzene.xyz > Benzindole.xyz
</pre><img src="Fusing.png"><pre>

Использование в режиме замещения:
$prog -ijk=2,11 -IJK=8,2 Indole.xyz Benzene.xyz > Phenylindole.xyz
</pre><img src="Substitution.png"><pre>
В режиме замещения номера k и K можно не указывать, и если они не указаны
(по крайней мере, K в заместителе), то будет действовать опция -rot. 

-R=substitutent  $prog -ijk=i,j,k -R=Me core.xyz
   Из __DATA__ зачитывается заместитель и помещается на место атома i.
   Расстояние между атомом j и первым атомом заместителя - сумма их 
     ковалентных радиусов.
-rot=degrees   Если задано, то заместитель поворачивается на degrees градусов.
   Одновалентность замещаемого атома не проверяется.
   
   Вместо -ijk можно задать -ii (только замещаемый атом), тогда j и k будут 
     определены автоматически. Дополнительные возможности -ii:
     -ii=m,n,... означает, что будут замещены атомы номер m,n,...
     -ii=element - будут замещены все элементы. Например, перфторировать:
     $0 -ii=H -R=F cores.xyz > coresF.xyz
 
-R (без указания заместителя) - будет напечатан список заместителей, 
    определенных в __DATA__, и сразу после этого выход.
    
-make_R=name  Сделать из молекулы заместитель. Нужно задать -ijk=i,j,k
   i - одновалентный атом, который потом будет удален из молекулы
   j - замещающий атом, к которому присоединен i
   k - еще какой-нибудь атом, задающий плоскость xy
   Например, $prog -make_R=Ph -ijk=7,1,2 Benzene.xyz
   Результат сохранить, записав его в __DATA__ (в конец файла программы $prog ).

HELP
  exit(0);
}

$rm_dist ||= 0.8;

if ($ii && $ijk) {
  die "Несовместимые опции -i и -ijk\n";
}
if ($ii && !$R) {
  die "-ii д.б. только для -R\n";
}

my ($i,$j,$k) = (1,2,3);
my $by_bond = 0;
if ($ijk) {
  $ijk=~/(\d+)[\s[:punct:]]+(\d+)(?:[\s[:punct:]]+(\d+))?/ ?
  ($i,$j,$k) = ($1,$2,$3) :
  die "Номера после -ijk должны быть целыми числами!\n";
  do {$by_bond = 1; $k = 0} unless $k;
  #die "Номера после -ijk не должны быть нулевыми!\n" if $i==0 or $j==0 or $k==0;
  die "Номера после -ijk не должны быть одинаковыми!\n" if $i==$j or $j==$k or $k==$i;
}
my ($I,$J,$K) = (1,2,3);
if ($IJK) {
  $IJK=~/(\d+)[\s[:punct:]]+(\d+)(?:[\s[:punct:]]+(\d+))?/ ?
  ($I,$J,$K) = ($1,$2,$3) :
  die "Номера после -IJK должны быть целыми числами!\n";
  $K ||= 0;
  #die "Номера после -IJK не должны быть нулевыми!\n" if $I==0 or $J==0 or $K==0;
  die "Номера после -IJK не должны быть одинаковыми!\n" if $I==$J or $J==$K or $K==$I;
}

if ($R && $R eq '1') {
  print 'Defined substitutents: ', join(' ', map {$_->[0]{Title}} &read_subst), "\n";
  exit;
}

if ($make_R) {
  my ($mol) = (read_molden())[-1];
  if ($i>$#$mol or $j>$#$mol or $k>$#$mol) {
    warn "Номер после -ijk больше числа атомов!\n";
    next;
  }
  if (colinearity($mol,$i,$j,$k)) {
    die "Atoms $i,$j,$k in the core lie on a single line\n";
  }
  re_orientation($mol,$i,$j,$k);
  splice @$mol,$i,1; # Удаляем Hi
  my $aj = splice @$mol, $j<$i ? $j : $j-1, 1; # Временно удаляем замещающий атом
  splice @$mol,1,0,$aj; # Вставляем его в начало
  $mol->[0]{Name} = $make_R;
  for (my $i=1; $i<@$mol; $i++) {
    $mol->[$i][1] *= -1;
  }
  write_molden($mol);
  exit;
}

my @mols = read_molden();
die "Замещение в 1- и 2-атомных молекулах не предусмотрено!\n" if grep {@$_<4} @mols;

# Ковалентные радиусы из Dalton Trans., 2008, 2832-2838
my %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  
);
  
if ($R) {
  my $degrees = defined($rot) ? $rot : '';
#  if ($R =~ s/,(-?\d+(?:\.\d+)?)$//) {
#    $degrees = $1;
#  }
  (my $subst0) = grep {$_->[0]{Title} eq $R} &read_subst;
  die "No $R in __DATA__\n" unless $subst0;
  # Будем немного поворачивать метилы, чтобы убрать симметрию
  if ($degrees eq '' && $R eq 'Me') {
    $degrees = 15;
  }
  $degrees ||= 0;
  # Зануляем координаты 1-го атома заместителя
  my ($x0,$y0,$z0) = ($subst0->[1][1],$subst0->[1][2],$subst0->[1][3]);
  for (my $i=1; $i<@$subst0; $i++) {
    $subst0->[$i][1] -= $x0;
    $subst0->[$i][2] -= $y0;
    $subst0->[$i][3] -= $z0;
  }
  # Поворачиваем
  if ($degrees) {
    rot_x($subst0,$degrees);
  }
  
  if ($ii) {
    # Парсим $ii
    if ($ii =~ /^[a-z][a-z]?$/i) {
      my $el = ucfirst lc $ii;
      $ii = '';
      for (my $n=1; $n<@{$mols[0]}; $n++) {
        $ii .= "$n," if $mols[0][$n][0] eq $el;
      }
      $ii =~ s/,$//;
      die "Нет элемента $el в 1-й структуре\n" unless $ii;
    }
    $ii =~ s/\.\./-/g;
    die "Invalid -ii=$ii\n" if $ii !~ /^\d+([,-]\d+)*/;
    my @ii;
    foreach (split /,/, $ii) {
      if (/(\d+)-(\d+)/) { # Раскрываем интервалы (i-n)
		    my @a = $1<=$2 ? $1..$2 : reverse $2..$1;
		    push @ii, @a;
	    }
      else {push @ii, $_}
    } 
    #print "@ii\n"; exit;
    
    foreach my $ii (@ii) {
      # Определяем соседей, т.к. они не заданы явно
      $i = $ii; # Замещаемый атом
      # Его ближайший сосед (в 1-й молекуле)
      $j = nearest($i, $mols[0]);
      # Ближайший сосед соседа
      $k = nearest($j, $mols[0], $i);
      
      # Расстояние между замещаемым атомом (i) и соседним к нему (j) в 1-й молекуле
      # Полагаем, что в остальных молекулах эти атомы такие же
      my $rij = dist($mols[0],$i,$j);
      
      # Локальная версия заместителя (будет модифицироваться)
      # $subst = $subst0; не годится (ссылки)
      my $subst = copy_mol($subst0); 
      #write_molden($subst);
      
      # Расстояние между соседним к замещаемому атомом и первым заместителя
      my $r = $radius{$mols[0][$j][0]}+$radius{$subst->[1][0]};
      # Сдвигаем заместитель по x, 
      # чтобы между соседним к замещаемому атомом и первым заместителя было $r
      for (my $i=1; $i<@$subst; $i++) {
        $subst->[$i][1] = $subst->[$i][1]+$rij-$r;
      }
      # Удаляем из заместителя 1-ю (служебную) строку
      shift @$subst;
      # 1-й атом заместителя
      my $first_atom = shift @$subst;
      my @mols_tmp;
      foreach  (@mols) {
        my $mol = copy_mol($_);
        # Подстраиваем молекулу под заместитель
        re_orientation($mol,$i,$j,$k);
        # Удаляем из молекулы замещаемый атом и на его место - 1-й атом заместителя
        splice @$mol, $i, 1, $first_atom;
        # Соединяем молекулу и заместитель без 1-го атома
        push @$mol, @$subst;
        push @mols_tmp, $mol;
      }
      @mols = @mols_tmp;
      #write_molden(@mols);
    }
  }
}
else {
  die "Нужно минимум две структуры\n" if @mols<2;

  my $subst = pop(@mols);
  $K = nearest($I,$subst,$J) if $K==0;
  #write_molden($subst); exit;
  if ($I>$#$subst or $J>$#$subst or $K>$#$subst) {
    die "Номер после -IJK больше числа атомов!\n";
  }
  if (colinearity($subst,$I,$J,$K)) {
    die "Atoms $I,$J,$K in the substitute lie on a single line\n";
  }
  re_orientation($subst,$I,$J,$K);
  if ($by_bond && $rot) {
    rot_x($subst,$rot);
  }

  foreach my $mol (@mols) {
    $k = nearest($i,$mol,$j) if $k==0;
    if ($i>$#$mol or $j>$#$mol or $k>$#$mol) {
      warn "Номер после -ijk больше числа атомов!\n";
      next;
    }
    if (colinearity($mol,$i,$j,$k)) {
      warn "Atoms $i,$j,$k in the core lie on a single line\n";
      next;
    }
    re_orientation($mol,$i,$j,$k);
    
    my @dist = distance($mol,$subst);
  #  # For debug
  #  for (my $i=1; $i<@$mol; $i++) {
  #    my @tmp;
  #    for (my $j=1; $j<@$subst; $j++) {
  #      push @tmp, sprintf("%d,%d=%.2f",$i,$j,$dist[$i][$j]) if $dist[$i][$j]<0.7;
  #    }
  #    print "@tmp\n" if @tmp;
  #  }
    
    for (my $i=1; $i<@$mol; $i++) {
      for (my $j=1; $j<@$subst; $j++) {
        if ($dist[$i][$j] < $rm_dist) {
          
          # Помечаем строчным шрифтом накладывающиеся водороды и др. одновалентные
          if ($mol->[$i][0] =~ /^(H|F|Cl|Br)$/ && $subst->[$j][0] !~ /^(H|F|Cl|Br)$/) {
            $mol->[$i][0] = lc($mol->[$i][0])
          }
          elsif ($mol->[$i][0] !~ /^(H|F|Cl|Br)$/ && $subst->[$j][0] =~ /^(H|F|Cl|Br)$/) {
            $subst->[$j][0] = lc($subst->[$j][0])
          }
          else {
          # Помечаем пробелом другие накладывающиеся атомы, в остове или в заместителе
          $rm_in_core ? $mol->[$i][0] = ' ' . $mol->[$i][0] :
                        $subst->[$j][0] = ' ' . $subst->[$j][0];
          }
        }
      }
    }
    if ($by_bond) {
      (my $at_m = $mol->[$i][0]) =~ s/\s+//g; $at_m = ucfirst lc $at_m;
      (my $at_s = $subst->[$J][0]) =~ s/\s+//g; $at_s = ucfirst lc $at_s;
      my $dx = $radius{$at_m}+$radius{$at_s}-$subst->[$J][1];
      for (my $i=1; $i<@$subst; $i++) {
        $subst->[$i][1] += $dx;
      }
    }
    push @$mol, @{$subst}[1..@$subst-1];
    $mol->[0] = '';
    @$mol = grep {$_ eq '' or $_->[0] !~ /^( |[a-z])/} @$mol;
    $mol->[0] = {};
  }
}
#pp @mols;
write_molden(@mols);


# Читает xyz. Параметры - имена xyz-файлов. Если параметров нет, то <>.
# Возвращает массив найденных молекул.
sub read_molden {
  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 = <>;
      if ($line =~ /^\s*(.*?)\s*$/) {
        $mol[0]{Title} = $1;
      }
      if ($line =~ /\s($num)\s/o) {
        $mol[0]{Energy} = $1;
      }
      if ($line =~ /Symmetry\s+(\w+)/) {
        $mol[0]{Symmetry} = $1;
      }
      #($mol[0]) = $line =~ /($num)/o;
      for (my $i=1; $i<=$N; $i++) {
        last LOOP if eof();
        next LOOP if eof(ARGV);
        $line = <>;
        #print $line;
        if ($line =~ /^\s*([A-Za-z]{1,2})\s+($num)\s+($num)\s+($num)\s*(.*)/io) {
          $mol[$i] = [$1,$2,$3,$4,$5];
          $mol[$i][0] = ucfirst(lc $mol[$i][0]);
        } else {
          next LOOP;
        }
      }
      push @mols, \@mol;
      last LOOP if eof();
    } else {
      undef $line;
    }
  }
  return @mols;
}

# Печатает xyz. Параметры -- список молекул, в конце м.б. имя файла.
# Если последний элемент списка - имя файла (не ссылка на массив),
# то печать в этот файл, иначе - на stdout
sub write_molden {
	my $fh = \*STDOUT;
  if (ref($_[-1]) ne 'ARRAY') {
    my $file = pop @_;
    open $fh, '>', $file or die "Can't write to $file: $!\n";
  }
  
  foreach my $mol (@_) {
		my $N = $#{$mol};
    print $fh " $N\n";
		print $fh "$mol->[0]{Name} " if exists $mol->[0]{Name};
		#print $fh " Symmetry $mol->[0]{Symmetry} " if $mol->[0]{Symmetry};
    print $fh "\n";
		for (my $i=1; $i<=$N; $i++) {
      my ($atom,$x,$y,$z) = @{$mol->[$i]};
      printf $fh " %-2s %12.8f %12.8f %12.8f", $atom, $x, $y, $z;
      #printf $fh uc($atom) eq 'H' ? " %10.3f" : " %9.2f"  , $ppm if $ppm;
      print $fh "\n";
		}
	}
  #close $fh;
}

# my ($x,$y,$z) = get_xyz($mol);
# Возвращает список ссылок на массивы 1..$N (@x,@y,@z) молекулы
# (нулевые элементы пустые)
sub get_xyz {
  my $mol = shift;
  my (@x,@y,@z);
  my $N = $#{$mol};
  for (my $i=1; $i<=$N; $i++) {
    $x[$i] = $mol->[$i][1];
    $y[$i] = $mol->[$i][2];
    $z[$i] = $mol->[$i][3];
  }
  return (\@x,\@y,\@z);
}

# put_xyz($mol,\@x,\@y,\@z);
# Помещает координаты из массивов 1..$N (@x,@y,@z) в молекулу
sub put_xyz {
  my ($mol,$x,$y,$z) = @_;
  my $N = $#{$mol};
  return undef if $#{$x} != $N || $#{$y} != $N || $#{$z} != $N;
  for (my $i=1; $i<=$N; $i++) {
    $mol->[$i][1] = $x->[$i];
    $mol->[$i][2] = $y->[$i];
    $mol->[$i][3] = $z->[$i];
  }
  1
}

# centre_to_atom($mol,$i)
# Помещает центр координат молекулы на атом $i
sub centre_to_atom {
  my ($mol,$i) = @_;
	my $N = $#{$mol};
  return undef if $i > $N;
	my ($xi,$yi,$zi) = ($mol->[$i][1],$mol->[$i][2],$mol->[$i][3]);
	for (my $n=1; $n<=$N; $n++) {
	  $mol->[$n][1] -= $xi;
    $mol->[$n][2] -= $yi;
    $mol->[$n][3] -= $zi; 
	}
  1
}

# ($A,$B,$C) = plane_normal($mol,$i,$j,$k);
# Возвращает вектор нормаль к плоскости, проходящей через атомы $i, $j, $k.
# Если атомов больше трех, плоскость через них проводится по наименьшим квадратам.
sub plane_normal {
	my ($mol,@nums) = @_;
  my ($x,$y,$z) = get_xyz($mol);
  my @x = @$x; my @y = @$y; my @z = @$z;
  my ($A,$B,$C);
	return undef if @nums<3;
	if (@nums == 3) {
	  my ($i,$j,$k) = @nums;
		#print "@nums\n";
    $A = ($y[$j]-$y[$i])*($z[$k]-$z[$i])-($z[$j]-$z[$i])*($y[$k]-$y[$i]);
		$B = ($z[$j]-$z[$i])*($x[$k]-$x[$i])-($x[$j]-$x[$i])*($z[$k]-$z[$i]);
		$C = ($x[$j]-$x[$i])*($y[$k]-$y[$i])-($y[$j]-$y[$i])*($x[$k]-$x[$i]);
		#print "@nums  $A $B $C\n";
	} else {
		my ($X,$Y,$Z,$XX,$YY,$ZZ,$XY,$YZ,$XZ);
		foreach (@nums) {
			$X  += $x[$_];        $Y  += $y[$_];        $Z  += $z[$_];
			$XX += $x[$_]*$x[$_]; $YY += $y[$_]*$y[$_]; $ZZ += $z[$_]*$z[$_];
			$XY += $x[$_]*$y[$_]; $YZ += $y[$_]*$z[$_]; $XZ += $z[$_]*$x[$_];
		}
		# A*XX + B*XY + C*XZ = -X   XX XY XZ | -X 	 -X XY XZ 	XX -X XZ	 XX XY -X
		# A*XY + B*YY + C*YZ = -Y   XY YY YZ | -Y 	 -Y YY YZ 	XY -Y YZ	 XY YY -Y
		# A*XZ + B*YZ + C*ZZ = -Z   XZ YZ ZZ | -Z 	 -Z YZ ZZ 	XZ -Z ZZ	 XZ YZ -Z
		$A = $X*($YZ*$YZ-$YY*$ZZ) + $Y*($XY*$ZZ-$YZ*$XZ) + $Z*($YY*$XZ-$XY*$YZ);
		$B = $X*($XY*$ZZ-$XZ*$YZ) + $Y*($XZ*$XZ-$XX*$ZZ) + $Z*($XX*$YZ-$XY*$XZ);
		$C = $X*($YY*$XZ-$YZ*$XY) + $Y*($YZ*$XX-$XY*$XZ) + $Z*($XY*$XY-$YY*$XX);
	}
  my $normal = sqrt($A**2 + $B**2 +$C**2);
  $A /= $normal; $B /= $normal; $C /= $normal;
	#print "@nums  $A $B $C\n";
	return ($A,$B,$C);
}

# re_orientation($mol,$i,$j,$k);
# Переориентирует молекулу так чтобы у атомов $i,$j,$k были координаты
# $i:  0.00  0.00  0.00
# $j:  f.ff  0.00  0.00
# $k:  f.ff  f.ff  0.00
# Dependecies: plane_normal, colinearity, put_xyz
sub re_orientation {
  my ($mol,$i,$j,$k) = @_;
	centre_to_atom($mol,$i);
	#&write_molden;
  my ($x,$y,$z) = get_xyz($mol);
  my @x = @$x; my @y = @$y; my @z = @$z;
 	
  my $xj = sqrt($x[$j]**2 + $y[$j]**2 + $z[$j]**2);
	my $xk = ($x[$k]*$x[$j] + $y[$k]*$y[$j] + $z[$k]*$z[$j]) / $xj;
	my $yk = sqrt($x[$k]**2 + $y[$k]**2 + $z[$k]**2 - $xk**2);
	my ($A1,$B1,$C1) = plane_normal($mol,$i,$j,$k);
  
	for (my $n=1; $n<=$#{$mol}; $n++) {
	  if ($n != $i && $n != $j && $n != $k) {
			my ($xn,$yn,$zn);
      if (colinearity($mol,$j,$k,$n)) {     # Try colinearity
        my $slope;
        my $xkj = $x[$k]-$x[$j];
        my $ykj = $y[$k]-$y[$j];
        my $zkj = $z[$k]-$z[$j];
        if (abs($xkj)>=abs($ykj) && abs($xkj)>=abs($zkj)) {
          $slope = ($x[$n]-$x[$j])/$xkj;
        }
        elsif (abs($ykj)>=abs($zkj)) {
          $slope = ($y[$n]-$y[$j])/$ykj;
        }
        else {
          $slope = ($z[$n]-$z[$j])/$zkj;
        }
        $zn = 0;
        $yn = $slope*$yk;
        $xn = $slope*($xk-$xj)+$xj;
      }
      else {
        $xn = ($x[$n]*$x[$j] + $y[$n]*$y[$j] + $z[$n]*$z[$j])/$xj;
        $yn = ($x[$n]*$x[$k] + $y[$n]*$y[$k] + $z[$n]*$z[$k] - $xn*$xk) / $yk;
        my $sqr = $x[$n]**2 + $y[$n]**2 + $z[$n]**2 - $xn**2 - $yn**2;
        $sqr = 0 if $sqr < 0;
        $zn = sqrt($sqr);
        my ($A2,$B2,$C2) = plane_normal($mol,$j,$k,$n);
        my $Vec = ($B1*$C2-$C1*$B2)*($x[$k]-$x[$j])+
                  ($C1*$A2-$A1*$C2)*($y[$k]-$y[$j])+
                  ($A1*$B2-$B1*$A2)*($z[$k]-$z[$j]);
        $zn = -$zn if $Vec < 0;
        #print "$Vec\n";
      }
			($x[$n],$y[$n],$z[$n]) =	($xn,$yn,$zn);
		}
	}
	($x[$j],$y[$j],$z[$j],$x[$k],$y[$k],$z[$k]) =	($xj,0,0,$xk,$yk,0);
	put_xyz($mol,\@x,\@y,\@z);
  #&write_molden;
}

sub colinearity {
  my ($mol,$i,$j,$k) = @_;
  my $eps = 1e-6;
  my $rij = sqrt(($mol->[$i][1]-$mol->[$j][1])**2+
                 ($mol->[$i][2]-$mol->[$j][2])**2+
                 ($mol->[$i][3]-$mol->[$j][3])**2);
  my $rik = sqrt(($mol->[$i][1]-$mol->[$k][1])**2+
                 ($mol->[$i][2]-$mol->[$k][2])**2+
                 ($mol->[$i][3]-$mol->[$k][3])**2);
  my $rjk = sqrt(($mol->[$k][1]-$mol->[$j][1])**2+
                 ($mol->[$k][2]-$mol->[$j][2])**2+
                 ($mol->[$k][3]-$mol->[$j][3])**2);
  my ($r1,$r2,$r3) = sort {$a<=>$b} ($rij,$rik,$rjk);
  if ($r1<$eps) {
    warn "Two atoms from $i,$j,$k coincide\n";
    return 1;
  }
  if ($r1+$r2-$r3<$eps) {
    return 1;
  }
  return undef;
}

sub distance {
  my ($mol,$subst) = @_;
  my @dist;
  for (my $i=1; $i<@$mol; $i++) {
    my ($x,$y,$z) = @{$mol->[$i]}[1..3];
    for (my $j=1; $j<@$subst; $j++) {
      my ($X,$Y,$Z) = @{$subst->[$j]}[1..3];
      $dist[$i][$j] = sqrt(($x-$X)**2+($y-$Y)**2+($z-$Z)**2);
    }
  }
  return @dist;
}

sub dist {
  my ($mol,$i,$j) = @_;
  my $rij = sqrt(($mol->[$i][1]-$mol->[$j][1])**2 +
                 ($mol->[$i][2]-$mol->[$j][2])**2 +
                 ($mol->[$i][3]-$mol->[$j][3])**2 );
  return $rij;
}

sub nearest {
 # Ближайший сосед к атому i (за исключением j, если он задан)
  my ($i,$mol,$j) = @_;
  my $N = @$mol;
  die "Index $i more then number of atoms ($N)\n" if $i > $N;
  my $k;
  my $r_min = 1e6;
   for (my $n=1; $n<$N; $n++) {
     next if $n==$i;
     next if $j and $n==$j;
     my $r = dist($mol,$i,$n);
     if ($r < $r_min) {
       $r_min = $r;
       $k = $n;
     }
   }
   return $k;
}

sub rot_x {
  my ($mol,$degrees) = @_;
  my $rad = $degrees/180*3.14159265358979;
  my $cos = cos($rad);
  my $sin = sin($rad);
  for (my $i=1; $i<@$mol; $i++) {
    my ($y,$z) = ($mol->[$i][2],$mol->[$i][3]);
    $mol->[$i][2] = $sin*$z+$cos*$y;
    $mol->[$i][3] = $cos*$z-$sin*$y;
  }
}

sub read_subst {
  my $tmp = "subst_DATA.tmp";
  open L, '>', $tmp or die "Can't write to file $tmp: $!\n";
  while (<DATA>) {
    print L;
  }
  close L;
  my @mols = read_molden($tmp);
  unlink $tmp;
  return @mols;
}

sub copy_mol {
  my $mol = shift;
  my $new_mol;
  my $N = $#{$mol};
  $new_mol->[0] = $mol->[0] if $mol->[0];
  for (my $i=1; $i<=$N; $i++) {
    $new_mol->[$i] = [@{$mol->[$i]}];
  }
  return $new_mol;
}

# Многоатомные заместители (-R) д.б. направлены вдоль связи в минус по x
# и замещающий атом д.б. первым
# Например, берем карбоксил C(1)O(2)OH. Делаем муравьиную кислоту H(5)C(1)O(2)OH, 
# reorient -ijk=1,5,2 и удаляем последний H
__DATA__ 
1
H
 H    0.00000000   0.00000000   0.00000000
1
F
 F    0.00000000   0.00000000   0.00000000
1
Cl
 Cl    0.00000000   0.00000000   0.00000000
1
Br
 Br    0.00000000   0.00000000   0.00000000
1
I
 I    0.00000000   0.00000000   0.00000000
2
OH
 O    0.00000000   0.00000000   0.00000000
 H   -0.36578731   0.89924971  -0.00160080
3
NH2
 N    0.00000000   0.00000000   0.00000000
 H   -0.46317533   0.84528872   0.31152275
 H   -0.46328438  -0.84608500   0.30949621
3
NO2
 N    0.00000000   0.00000000   0.00000000
 O   -0.58263507  -1.09524818   0.00003862
 O   -0.58242824   1.09527120  -0.00021915
6
Ac
 C    0.00000000   0.00000000   0.00000000
 O   -0.63621472   1.05590662  -0.01044610
 C   -0.72946410  -1.32068058   0.01263019
 H   -1.80838042  -1.13662904   0.01920459
 H   -0.47664668  -1.88329053   0.91498375
 H   -0.48874032  -1.89425334  -0.88611337
7
OAc
 O    0.00000000   0.00000000   0.00000000
 C   -2.09330995  -0.51410490   0.90453115
 C   -0.60525311  -0.59682033   1.07615371
 O   -0.06128596  -1.12023185   2.03931056
 H   -2.58431301  -0.97683133   1.76547623
 H   -2.39259641  -1.05232417   0.00170809
 H   -2.40264923   0.53269479   0.85010802
7
Tf
 S    0.00000000   0.00000000   0.00000000
 O   -0.46990098   1.25129727   0.58378706
 O   -0.46963891  -1.25449179   0.57725421
 C   -0.32886683   0.00453706  -1.78988747
 F    0.19848152   1.09993472  -2.40786790
 F    0.20192515  -1.08562582  -2.41423083
 F   -1.66239701   0.00318499  -2.06980455
7
OAcf
 O    0.00000000   0.00000000   0.00000000
 C   -0.56097445   0.02867405  -1.26068702
 O    0.05438301   0.05268514  -2.32316466
 C   -2.07555897   0.02563674  -1.11649931
 F   -2.53069833   1.10468748  -0.41666554
 F   -2.53061486  -1.07491367  -0.45104008
 F   -2.73184914   0.04468208  -2.31039581
8
OTf
 O    0.00000000   0.00000000   0.00000000
 S   -0.93299857   1.34153210   0.00648950
 O   -0.72457998   2.03808068   1.26715589
 O   -0.73477125   2.04327390  -1.25297271
 C   -2.54432913   0.48967807   0.01138638
 F   -3.59061694   1.36309622   0.03407336
 F   -2.71519003  -0.29706564  -1.09003834
 F   -2.69303641  -0.32724073   1.09398841
5
OSO2F
 O    0.00000000   0.00000000   0.00000000
 S   -0.92591248   1.33942816  -0.00447774
 O   -0.76864144   2.02962357  -1.26910829
 O   -0.78109517   2.02890033   1.26210959
 F   -2.43960961   0.79721767  -0.01189763
4
Me
 C    0.00000000   0.00000000   0.00000000
 H   -0.36633733   1.03626161   0.00000000
 H   -0.36654080  -0.51811160  -0.89737613
 H   -0.36644444  -0.51815003   0.89737614
7
Et
 C    0.00000000   0.00000000   0.00000000
 H   -0.36504683  -0.89498626   0.53324939
 H   -0.36497446   0.87692574   0.56249613
 C   -0.58915384   0.02348760  -1.42230384
 H   -1.69099508   0.02433998  -1.39139030
 H   -0.25883522   0.92167552  -1.96832362
 H   -0.26095308  -0.85749903  -1.99687644
10
i-Pr
 C    0.00000000   0.00000000   0.00000000
 C   -0.57117591   0.65646117  -1.26095149
 C   -0.57440975   0.65848134   1.25838770
 H   -0.35398868  -1.04062877   0.00063689
 H   -1.66134362   0.55232549  -1.28627761
 H   -0.34005114   1.72626932  -1.30559252
 H   -0.17014059   0.18423811  -2.16451648
 H   -1.66478093   0.55345206   1.28160660
 H   -0.34454695   1.72857276   1.30193579
 H   -0.17565424   0.18831549   2.16392337
13
t-Bu
 C    0.00000000   0.00000000   0.00000000
 C   -0.48577368  -0.12183199   1.45840990
 C   -0.60894131   1.28909972  -0.60458711
 C   -0.59218571  -1.17598353  -0.81566374
 H   -1.58061646  -0.12263987   1.51246765
 H   -0.12570329   0.71400367   2.06971604
 H   -0.13267812  -1.04943249   1.92370931
 H   -1.70126324   1.21732095  -0.67157203
 H   -0.22927260   1.47094779  -1.61659307
 H   -0.39700901   2.17244297   0.00741868
 H   -1.68403037  -1.09971749  -0.88483855
 H   -0.19694856  -1.18725041  -1.83796447
 H   -0.38301810  -2.14807743  -0.35655940
10
Pr
 C    0.00000000   0.00000000   0.00000000
 C   -0.57253822   0.01328753   1.41944046
 C   -2.09299017  -0.00217992   1.40571899
 H   -0.37170184   0.86886507  -0.55928362
 H   -0.37232757  -0.87887159  -0.54313761
 H   -0.21137897  -0.85722960   1.98004913
 H   -0.22970179   0.90623024   1.95584359
 H   -2.48991413   0.87520273   0.88530192
 H   -2.48071472   0.00587743   2.42910044
 H   -2.47167279  -0.89985000   0.90662268
13
Bu
 C    0.00000000   0.00000000   0.00000000
 C   -0.57360691  -0.01947928   1.42012259
 C   -2.10241636  -0.00781196   1.40006447
 C   -2.67808513  -0.03021399   2.80707062
 H   -0.37162490   0.88155696  -0.53899430
 H   -0.37217276  -0.86629993  -0.56278915
 H   -0.22182412  -0.91287419   1.95150042
 H   -0.20953414   0.85080336   1.98064781
 H   -2.47494048  -0.87655317   0.84504174
 H   -2.46032447   0.88810527   0.87979525
 H   -3.77174006  -0.00725817   2.76942039
 H   -2.33946895   0.83646218   3.38380497
 H   -2.37539023  -0.93729828   3.33986305
5
CH2OH
 C    0.00000000   0.00000000   0.00000000
 O   -0.49630939  -0.70500320  -1.13264246
 H   -0.38639904  -0.48020672   0.90523548
 H   -0.40889541   1.01524443  -0.03595079
 H    0.07185703  -0.46122780  -1.88399755
5
OMe
 O    0.00000000   0.00000000   0.00000000
 C   -0.66956608  -1.25374046   0.00960609
 H   -1.74622758  -1.05970707   0.01455469
 H   -0.42659423  -1.82134099   0.91379745
 H   -0.43651836  -1.82978356  -0.89174250
4
CF3
 C    0.00000000   0.00000000   0.00000000
 F   -0.52144780  -1.08032120  -0.63991068
 F   -0.52063277   1.09753291  -0.61070772
 F   -0.49918527  -0.01671036   1.26320528
9
NMe2
 N    0.00000000   0.00000000   0.00000000
 C   -0.61774579  -1.21677782   0.52254975
 C   -0.61727376   1.19677460   0.56753717
 H   -1.71141762  -1.14876426   0.48744881
 H   -0.32262348  -1.40629297   1.56031517
 H   -0.36208023  -2.08826095  -0.08796591
 H   -1.71102230   1.13008490   0.53095136
 H   -0.32121346   1.34786586   1.61138665
 H   -0.36214312   2.09032342  -0.01047449
6
NHMe
 N    0.00000000   0.00000000   0.00000000
 C   -0.70019673   1.11953005   0.59824877
 H   -0.39652901  -0.89636421   0.27035688
 H   -1.77897029   0.93536684   0.56074501
 H   -0.52270908   2.05218844   0.05387282
 H   -0.42268500   1.25218243   1.64906909
4
SO2F
 S    0.00000000   0.00000000   0.00000000
 O   -0.49967280  -1.27599923   0.48263179
 O   -0.49972807   1.28787352   0.44986049
 F   -0.24464621  -0.02088426  -1.62430688
5
SO2H2
 S    0.00000000   0.00000000   0.00000000
 O   -0.19595455   1.16046815   1.13338642
 O   -0.21196581  -1.44107391   0.74563706
 H   -1.08932176   1.57173901   1.02582623
 H   -1.11580778  -1.78055834   0.53241897
13
SiMe3
 Si   0.00000000   0.00000000   0.00000000
 C   -0.63147386  -1.74918443   0.00471649
 C   -0.54039508   0.89640044  -1.54115557
 C   -0.54052888   0.90477339   1.53664156
 H   -1.72592624  -1.76363375   0.00336355
 H   -0.29194381  -2.29707491  -0.87960166
 H   -0.29402492  -2.29193898   0.89299893
 H   -1.63212810   0.92071191  -1.61227560
 H   -0.18173646   1.93014167  -1.54578721
 H   -0.15649590   0.40061656  -2.43796010
 H   -1.63199732   0.91552354   1.61574437
 H   -0.14344544   0.42320469   2.43549151
 H   -0.19533199   1.94311668   1.52820068
4
SiCl3
 Si   0.00000000   0.00000000   0.00000000
 Cl  -0.61928575  -0.02601367  -1.93199234
 Cl  -0.70519761   1.66425977   0.92481108
 Cl  -0.69890157  -1.64352683   0.96658154
5
OCF3
 O    0.00000000   0.00000000   0.00000000
 C   -0.63172280  -0.17883524   1.23310657
 F   -1.93781571  -0.21795343   0.98423981
 F   -0.25051234  -1.32188671   1.82809634
 F   -0.36433675   0.83687844   2.07236479
25
Ad
 C    0.00000000   0.00000000   0.00000000
 C    1.54198334   0.00000000   0.00000000
 C   -0.51133498   1.45472796   0.00000000
 C   -0.51157620  -0.72201904   1.26278090
 C    2.06101588   0.72853793   1.25601024
 C    1.54566008   0.00513332   2.51644614
 C    0.00371049   0.00517810   2.52107383
 C    1.54583374   2.18190863   1.25366149
 C    0.00385296   2.18596390   1.25603170
 C   -0.50766549   1.45986852   2.51647342
 H    1.92001424  -1.03893166  -0.01771768
 H    1.91995436   0.50061616  -0.91055943
 H   -0.16511455  -1.77204384   1.26436657
 H   -1.61709255  -0.74294800   1.26441490
 H   -0.36394009  -0.51376145   3.42536585
 H    1.92643650   0.50928722   3.42388373
 H    1.92370460  -1.03373968   2.53718813
 H    3.16655244   0.72757383   1.25436565
 H    1.92367950   2.71530006   0.36185054
 H    1.92673074   2.71960815   2.14156998
 H   -0.36369185   3.22861167   1.25440295
 H   -1.61682934   1.46694509  -0.01785233
 H   -0.16416207   1.97732192  -0.91046922
 H   -0.15829637   1.98628418   3.42390724
 H   -1.61312020   1.47198756   2.53725220
9
HOTf1
 H    0.00000000   0.00000000   0.00000000
 O   -0.98162922   0.00000000   0.00000000
 S   -1.45258609   1.57384489   0.00000000
 O   -1.10467936   2.11530184  -1.29820541
 O   -1.10991783   2.17103437   1.28442734
 C   -3.33263188   1.30340313   0.00000108
 F   -3.70120578   0.63398927  -1.09655970
 F   -3.68043493   0.60277027   1.09239431
 F   -3.93292072   2.49932250   0.02480991
9
HOTf
 H   -0.00000000   0.00000000   0.00000000
 O   -0.66577160   0.72134880   0.00000050
 S   -2.14172689   0.00000000   0.00000000
 O   -2.30365453  -0.62289267   1.29820498
 O   -2.34816239  -0.65684095  -1.28442779
 C   -3.21809895   1.56497096   0.00000000
 F   -2.97615948   2.28983359   1.09656128
 F   -2.93913081   2.29574538  -1.09239273
 F   -4.50405320   1.19498248  -0.02480909
7
COOMe
 C    -1.08900009   0.00000000   0.00000000
 O    -1.61230432   1.09602388   0.00000000
 O    -1.76278306  -1.16702606   0.00000006
 C    -3.20233111  -1.02342070  -0.00539596
 H    -3.59269073  -2.04765642  -0.01221806
 H    -3.53979578  -0.48307191   0.89138050
 H    -3.53177852  -0.47348547  -0.89922268
11
Ph 
 C    -1.08900029   0.00000000   0.00000000
 C    -1.78900010   1.21243567   0.00000000
 C    -3.18900048   1.21243567   0.00000000
 C    -3.88900018   0.00000050   0.00000000
 C    -3.18900036  -1.21243518   0.00000000
 C    -1.78900085  -1.21243567   0.00000000
 H    -1.24449967   2.15553742   0.00000000
 H    -3.73350034   2.15553742   0.00000000
 H    -4.97800047   0.00000050   0.00000000
 H    -3.73350080  -2.15553693   0.00000000
 H    -1.24450100  -2.15553742   0.00000000