#!/usr/bin/perl -ws

#use Data::Dump 'pp';
our ($l,$n,$h,$help);

# Если опция -h, печатаем справку
if ($h || $help) {
  (my $program = $0) =~ s/^.*[\/\\]//;
  print "
Принадлежности колебания той или иной части молекулы.
Молекула делится на части и вычисляется, насколько колебание изменяет каждую 
часть. Численным критерием является величина с размеренностью момента инерции:
RMS = sum(M*(r-r0)**2), где M - атомная масса, (r-r0) - смещение атома в 
результате колебания (bohr), суммирование по всем атомам данной части молекулы.


Usage: $program [-l=min,max] part1_part2... filename
Зависимости: perl

filename - файл, в котором д.б. равновесная геометрия и колебания в виде
           смещений атомов по x,y,z.
           Автоматически распознаются форматы: freq-файл molden'а и пока все.
           Файлы в таком формате могут быть созданы квантово-химическим
           визуализатором molden (Write--Molden format), либо из расчета 
           Природой скриптом pri2mol.
part1 - номера атомов части молекулы через запятую, можно с интервалами 
        через дефис, например 2-5,7-9 то же самое что 2,3,4,5,7,8,9
        Части разделяются подчерком или пробелом (в последнем случае с 
        кавычками \"1-5 6-10\").
part2 - то же самое для второй части молекулы и т.д.

-l=min,max  Пределы RMS для статистики. По умолчанию -l=0.1,0.9
-n          Нормировать ли RMS частей на RMS целой молекулы

В результате выдается табличка со следующими колонками:
Vibr       - номер колебания
Freq       - частота колебания
Intens     - интенсивность колебания (если есть в файле)
Whole_vibr - величина RMS для целой молекулы
Part 1     - величина RMS для части 1
Part 2     - величина RMS для части 2
....

$program -help дает более подробное описание.
  ";
  if ($help) {
    print "
    Написание этой программы инициировано обсуждением на chemport'е:
       http://www.chemport.ru/forum/viewtopic.php?f=11&t=49893, 
       там есть пример с разъяснениями.
    
    Некоторые замечания.
       Алгоритм вычисления величины RMS такой: берется часть равновесной
    геометрии (атомов, заданных декартовыми координатами), перемещается в центр
    масс. На нее накладывается колебание  (смещения атомов по x,y,z, тоже на 
    всякий случай перенесенные в их \"центр масс\"). Новая геометрия максимально
    совмещается с равновесной  (критерий - минимальная RMS между ними). Эта RMS
    и будет характеризовать размах колебания атомов данной части.
       В квантово-химических программых (по крайней мере в GAMESS)
    колебания нормированы на 1 ед. величины RMS. Поэтому для целой молекулы  в
    случае \"нормальных\" колебаний RMS между равновесной и новой геометриями
    равно 1 (или близко к тому, хотя и не всегда). А для трансляций/вращений,
    молекула перемещается как целое, и RMS получается, естественно, 0. Если для
    одной  части RMS близко к 0, а для другой - к ~1, то значит, первая часть в 
    колебания перемещается как целое, а вторая искажается. В случае нулей для
    всех частей (и 1 для целой молекулы) фрагиенты колеблются друг относительно
    друга как целые.
       Молекулу не обязательно разбивать на две непересекающиеся части. Можно 
    указать только одну часть, или три и больше (например, указать атомы 
    карбонильной группы, чтобы найти ее характеристическое колебание (зачем 
    только это могло бы понадобиться? Разве что для автоматизации)).
       Можно давать в качестве параметров несколько файлов, но для всех них
    разбиение на части будет одно и то же.
    \n";
  }
  exit;
}

# Массы изотопов
my %massa = qw(
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
);
#my $BA = 0.529177249;

$ARGV[0] =~ s/_/ /g;
my @V = map {[get_atom_numbers($_)]} split ' ', shift;
#pp(\@v); exit;
for (my $i=0; $i<@V; $i++) {
  print "Part ", $i+1, ": atoms @{$V[$i]}\n";
}
print "\n";

my $min = 0.1;
my $max = 0.9;
if ($l && $l =~ /^(\d+\.\d+),(\d+\.\d+)$/) {
  ($min,$max) = ($1,$2);
}

foreach my $file (@ARGV) {
  my ($mol,$freq,$intens,$vibr,$outtype) = read_IR_file($file);
  if ($outtype ne 'molden') {
    print "Not Molden format of file $file\n";
    print "Try convert to Molden format via molden (Write button)\n";
    next;
  }
  print "File $file of $outtype format\n";
  
  foreach (@V) {
    if (grep {$_>$#{$mol}} @$_) {
      die "There is atom number > $#{$mol}. Die.\n";
    }
  }
  
  die "Freq != Vibr\n" if @$freq != @$vibr;
  for (@$vibr) {die "mol != vibr\n" if @$_ != @$mol};
  
  my @v = @V;
  unshift @v, [1..$#{$mol}];
  #pp(\@v); exit;
  
  my @mol0;
  foreach (@v) {
    my $mol0 = [{}, map {$mol->[$_]} @$_];
    $mol0 = copy_mol($mol0);
    to_centre_mass($mol0);
    push @mol0, $mol0;
    #pp($mol0);
  }
  
  #pp(\@mol0); exit;
  
  my @rms;
  for (my $i=0; $i<@$freq; $i++) {
    my @mol1;
    for (my $j=0; $j<@v; $j++) {
      my $mol1 = [{}, map {$vibr->[$i][$_]} @{$v[$j]}];
      $mol1 = copy_mol($mol1);
      to_centre_mass($mol1);
      for (my $k=1; $k<@$mol1; $k++) {
        $mol1->[$k][1] += $mol0[$j][$k][1];
        $mol1->[$k][2] += $mol0[$j][$k][2];
        $mol1->[$k][3] += $mol0[$j][$k][3];
      }
      push @mol1, $mol1;
    }
    #pp(\@mol1); exit;
    
    # массив RMS. В 0-м элементе макс. RMS
    $rms[0] = 0;
    for (my $j=0; $j<@v; $j++) {
      $rms[$i+1][$j] = get_rms($mol0[$j],$mol1[$j]);
      $rms[0] = $rms[$i+1][$j] if $rms[0] < $rms[$i+1][$j];
    }
    #printf "%10.2f\n", get_rms($mol,$vibr->[$i]);
    #pp @rms; 
  }
  
  # Нормировка
  if ($n) {
    for (my $i=1; $i<@rms; $i++) {
      my $f = $rms[$i][0];
      next if $f < 0.1;
      foreach (@{$rms[$i]}) {
        $_ /= $f;
      }
    }
  }
  
  # Печать
  printf "%4s%10s%10s    %10s", qw(Vibr Freq Intens Whole_vibr);
  for (my $i=0; $i<@V; $i++) {
    printf "%10s", "Part ".($i+1);
  }
  print "\n";
  
  for (my $i=1; $i<@rms; $i++) {
    printf "%4d%10.2f%10.2f    ", $i,$freq->[$i-1]||0,$intens->[$i-1]||0;
    foreach (@{$rms[$i]}) {
      printf "%10.2f", $_;
    }
    print "\n";
  }
  
  ## Statistics ##
  my @trans = grep {$rms[$_][0]<$min} 1..$#rms;
  print "\n";
  print "translations/rotayions of whole mol: " ,join(',',@trans), "\n" if (@trans);
  
  my %t = map {$_=>1} @trans;
  my @vibr = grep {! $t{$_}} 1..$#rms;
  
  print "Norm. RMS of the parts on whole mol. RMS\n";
  
  my (%nul,%one,$yes);
  for (my $j=1; $j<@{$rms[1]}; $j++) {
    foreach (@vibr) {
      my $rms = $rms[$_][$j];
      $rms /= $rms[$_][0] if $n;
      push @{$one{$_}},$j if $rms >= $max;
      push @{$nul{$_}},$j if $rms <= $min;
      push @{$yes{$_}},$j if $rms > $min && $rms < $max;
    }
  }
  #pp(\%nul, \%yes, \%one);
  
  $nul{max} = length "RMS<$min"; 
  $yes{max} = length "$min<RMS<$max"; 
  $one{max} = length "RMS>$max";
  foreach (@vibr) {
    if (exists $nul{$_}) {
      $nul{$_} = join ',', @{$nul{$_}};
      $nul{max} = length($nul{$_}) if $nul{max} < length($nul{$_});
    } else {$nul{$_} = ''}
    if (exists $yes{$_}) {
      $yes{$_} = join ',', @{$yes{$_}};
      $yes{max} = length($yes{$_}) if $yes{max} < length($yes{$_});
    } else {$yes{$_} = ''}
    if (exists $one{$_}) {
      $one{$_} = join ',', @{$one{$_}};
      $one{max} = length($one{$_}) if $one{max} < length($one{$_});
    } else {$one{$_} = ''}
  }
  
  
  printf "\n%-5s | %-$nul{max}s | %-$yes{max}s | %-$one{max}s\n", 
         'Vibr.', "RMS<$min", "$min<RMS<$max", "RMS>$max";
  print '-'x(14+$nul{max}+$yes{max}+$one{max}), "\n";
  foreach (@vibr) {
    printf "%-5d | %-$nul{max}s | %-$yes{max}s | %-$one{max}s\n",
    $_, $nul{$_}, $yes{$_}, $one{$_};
  }
}

sub get_atom_numbers {
	my ($vector,@result);
	$vector = shift;
	$vector =~ s/^\s+//;   # Убираем пробелы в начале
	$vector =~ s/\s+$//;   # Убираем пробелы в конце
	$vector =~ s/\s+/,/g;  # Заменяем пробелы на запятые

	unless ($vector =~ /^(\d+(-\d+)?,)*\d+(-\d+)?$/g) {
  	die "Not valid atom numbers: $vector\n";
	}

	foreach (split /,/, $vector) {
	  push @result, /(\d+)-(\d+)/ ? ($1<=$2 ? $1..$2 : reverse($2..$1)) : $_;
	}
	return sort {$a<=>$b} @result;
}

sub to_centre_mass {
# centre_to_atom($mol)
# Помещает центр координат молекулы в центр масс
  my $mol = shift;
	my $N = $#{$mol};
  #return undef if $i > $N;
	
  my ($x0,$y0,$z0,$M);
  for (my $i=1; $i<=$N; $i++) {
	  my $m = $massa{$mol->[$i][0]};
    $x0 += $m * $mol->[$i][1];
	  $y0 += $m * $mol->[$i][2];
	  $z0 += $m * $mol->[$i][3];
    $M += $m;
	}
  $x0 /= $M; 
  $y0 /= $M; 
  $z0 /= $M; 
  
	for (my $i=1; $i<=$N; $i++) {
	  $mol->[$i][1] -= $x0;
    $mol->[$i][2] -= $y0;
    $mol->[$i][3] -= $z0; 
	}
#	my $rms;
#  for (my $i=1; $i<=$N; $i++) {
#	  $rms += $massa{$mol->[$i][0]}*($mol->[$i][1]**2+$mol->[$i][2]**2+$mol->[$i][3]**2);
#	}
#  return sqrt($rms);
  return 1;
}

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-9;
  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 (abs($phix) < $FLT_EPSILON &&
	           abs($phiy) < $FLT_EPSILON &&
	           abs($phiz) < $FLT_EPSILON &&
	           ($e == $e0 || abs(($e-$e0)/($e+$e0)) < $FLT_EPSILON));
  }

  #/*** 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 sum {
  my $sum = 0;
  foreach (@_) {
    $sum += $_;
  }
  return $sum;
}

sub read_IR_file {
	my $file = shift;
  my $d = qr/\d+(?:\.\d+)?/;
  my $outtype = 'unknown';
  my $debug = 0;
  
  open PRI, $file or die "Can't open $file: $!\n";
  # determine the creator of the file
  while (<PRI>) {
    if (m/Priroda(?:\s+version)?\s+(\d+)/) { # for Priroda v. 2-9
      my $ver = $1;
      last if eof(PRI); $_ = <PRI>;
      next unless m/^\s*copyright\s+\(c\)\s.*\sLaikov$/; # for Priroda v. 2-9
      $outtype = "priroda$ver";
      last;
    }
    elsif (/^\s{8,}\*{30,}$/) {
      last if eof(PRI); $_ = <PRI>;
      next unless m/^\s*\*\s+GAMESS\s+VERSION\s+=\s+.+\s+\*$/;
      last if eof(PRI); $_ = <PRI>;
      next unless m/^\s*\*\s+FROM\s+IOWA\s+STATE\s+UNIVERSITY\s+\*$/;
      $outtype = 'gamessUS';
      last;
    }
    elsif (m/\Q[Molden Format]/) {
      $outtype = 'molden';
      last;
    }
    elsif (m/^\s*$d\s+$d\s*$/) {
      $outtype = 'txt';
      seek PRI, 0, 0;
      last;
    }
  }
  warn "$file: $outtype\n" if $debug;
  
  my (@mol,@freq,@intens,@vibr);
	
  ### Parse priroda 4,6 format
  if ($outtype eq 'priroda6') {
	  while (<PRI>) {
		  if (m/\Q | Mode | Freq.   | Mass. | IR Int. /) { # Priroda 6
		    <PRI>;
        while (<PRI>) {
		      last if m/\*/;
          s/\|//g;
          my ($freq,$intens) = (split)[1,3];
          $freq =~ s/i//;
          push @freq, $freq;
          push @intens, $intens;
		    }
  		}
	  }
	}
  
  ### Parse priroda2 format
  elsif ($outtype eq 'priroda2') {
    while (<PRI>) {
      push @freq, grep {!/^i/} split if s/^ freq\.//;
		  push @intens, split if s/^ ir i\.//;
    }
  }
  
  ### Parse gamess US format
  elsif ($outtype eq 'gamessUS') {
	  while (<PRI>) {
      push @freq, grep {!/^I/} split if s/^       FREQUENCY://;
		  push @intens, split if s/^    IR INTENSITY://;
	  }
  }
  
  ### Parse molden freq format
  elsif ($outtype eq 'molden') {
    my ($key,$i);

    while (<PRI>) {
      if (m/^\[(.*?)\]/) {
        $key = $1;
        next;
      }
      next unless defined $key;

      chomp;
      if ($key eq 'FREQ') {
        push @freq, 0+$_;
      }
      elsif ($key eq 'INT') {
        push @intens, 0+$_;
      }
      elsif ($key eq 'FR-COORD') {
        push @mol, [split ' '];
      }
      elsif ($key eq 'FR-NORM-COORD') {
        $i=$1,next if m/^vibration\s+(\d+)/;
        push @{$vibr[$i-1]}, [split ' '];
      }
    }

    
    unshift @mol, {};
    foreach (@vibr) {
      unshift @$_, {};
      for (my $i=1; $i<@mol; $i++) {
        unshift @{$_->[$i]}, $mol[$i][0];
      }
    }
  }
  
  elsif ($outtype eq 'txt') {
	  while (<PRI>) {
	    m/^\s*($d)\s+($d)\s*$/ or last;
      push @freq,   $1;
      push @intens, $2;
	  }
  }
    
	close PRI;
  
#  if ($#freq != $#intens) {
#		warn "freq != intens for $file";
#		return '';
#	}
  
  #pp(\@vibr); exit;
  if ($debug) {
    for (my $i=0; $i<@freq; $i++) {
      print "$freq[$i]\t$intens[$i]\n";
    }
  }
  
  return \@mol,\@freq,\@intens,\@vibr,$outtype;
}

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