Лаборатория изучения механизмов органических реакций

renum

#!/usr/bin/perl -ws

our($h,$help,$p,$prN,$crop,$canon);

# Если опция -h, печатаем справку
if ($h or $help) {
  ($program = $0) =~ s/^.*[\/\\]//;
  print "
Перенумерация xyz-файлов

Usage: $program vector [filename] [filenames] [> results]

Зависимости: perl, [Chemistry::Mol, Chemistry::Canonicalize, obabel].

filename - файл с декартовыми координатами в формате molden'а 
           (если отсутсвует, то ожидается ввод с STDIN).
           Генерится программами <pri2mol> или <g2i>
vector - новая последовательность атомов (вначале перечисленные,
         затем остальные по-порядку), например
         6,2-5 - поменять местами 1-й и 6-й атомы (6,2,3,4,5,1,7,...)
         '9,5' - эквивалентно 9,5,1,2,3,4,6,7,8,10,...
         3-1   - эквивалентно 3,2,1
Опции:
-crop : сгенеровать xyz-файл только из атомов, перечисленных в vector
        (в порядке, заданном vector'ом)
-p='7,8 3,4 6,9': пары атомов 7 и 8, 3 и 4, 6 и 9 поменять местами
                  vector с этой опцией не нужен
-auto : автоматическая перенумерация с помощью obabel. Образец для 
        перенумерации берется из 1-й молекулы, т.е. ее нумерация не меняется.
        При -auto=0 все молекулы будут иметь канонические номера 
        (весьма странные). vector с этой опцией не нужен.
-canon : автоматическая перенумерация с помощью Chemistry::Canonicalize.
         При -canon или -canon=1 Образец для перенумерации берется из 
         1-й молекулы, т.е. ее нумерация не меняется.
         При -canon=0 все молекулы будут иметь канонические номера 
         (весьма странные). vector с этой опцией не нужен.
-prN :  допечатывать в output номера атомов, которые были в input'е
\n";
  exit;
}

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

BEGIN {
  if (defined $canon) {
    unless (eval "use Chemistry::Mol; 1") {
      die "Coudn't load Chemistry::Mol: $@\n";
    }
    unless (eval "use Chemistry::Canonicalize 'canonicalize'; 1") {
      die "Coudn't load Chemistry::Canonicalize: $@\n";
    }
  }
}

my $obabel_x = '/usr/local/bin/obabel';
if (defined $auto) {
  unless (-x $obabel_x) {
    die "Can not run $obabel_x which needs for -auto option\n";
  }
  obabel_canonical(@ARGV);
  exit;
}

unless (defined($canon) or defined($auto)) {
  if ($p) {
    my @pairs = split ' ', $p;
    my @vector;
    foreach (@pairs) {
      m/^(\d+),(\d+)$/ || die "Invalid pair $_\n";
      ($vector[$2],$vector[$1]) = ($1,$2);
    }
    for (my $i=1; $i<@vector; $i++) {
      $vector[$i] = $i unless $vector[$i];
    }
    shift @vector;
    $vector = join ',', @vector;
  }
  else {
    $vector = shift;
  }

  $vector =~ s/^\s+//;      # Убираем пробелы в начале
  $vector =~ s/\s+$//;      # Убираем пробелы в конце
  $vector =~ s/[\s,]+/,/g;  # Заменяем пробелы на запятые

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

  # Раскрываем интервалы (i-n)
  #while ($vector =~ /((\d+)-(\d+))/) {
  #  my ($min,$max,$ttt);
  #  die "Error in permutation vector: $2-$3\n" if $2==$3;
  #  if ($2<$3) {
  #    for ($2..$3-1) {$ttt .= "$_,"};
  #  } else {
  #    for (my $i=$2; $i>$3; $i--) {
  #      $ttt .= "$i,";
  #    }
  #  }
  #   $ttt .= $3;
  #  $vector =~ s/$1/$ttt/;
  #}
  #
  #@vector = split /,/, $vector;

  foreach (split /,/, $vector) {
    if (/(\d+)-(\d+)/) { # Раскрываем интервалы (i-n)
      my @a = $1<=$2 ? $1..$2 : reverse $2..$1;
      push @vector, @a;
    }
    else {push @vector, $_}
  } 

  #print "@vector\n"; exit;

  # Находим максимальный номер
  $max = 0;
  foreach (@vector) {
    die "Error: atom number == 0\n" if $_<=0;
    warn "Warning! The same numbers: $_\n" if $vector{$_}++; # Проверка на одинаковые номера
    $max = $_ if $max < $_;
  }
  #if ($max != $#vector+1) {print "Error: max atom number is't true\n"; exit};

  unless ($crop) {
    # Дописываем массив номеров до максимального (5,8 => 5,8,1,2,3,4,6,7)
    foreach (1..$max) {
      undef @is_this_number;
      for (@vector) { $is_this_number[$_] = 1; }
      push @vector, $_ unless $is_this_number[$_];
    }
  }

  #foreach (@vector) {print "$_ "}; print "\n";
  #exit;
}

my $count;
#my (@sample_vector,@curr_vector); # for -canon
# Читаем input
while (<>) {
  if (/^\s*(\d+)\s*$/) {
    $N = $1;
    $crop ? print(' ', scalar(@vector), "\n") : print;
  }  else {
#    print "Not molden format of input";
    next;
  }

  # 2-nd string.
  $_ =<>;
  print;
  
  # Читаем геометрию с свойства
  for ($i=1;$i<=$N;$i++) {
    $atom[$i] = <>;
  }
  chomp @atom;
  
  die "Error: too big max atom number: $max\n" if (!defined($canon) && $max > $N);

  $count++;
  if (defined $canon) {
    if ($canon == 0) {
      @vector = canon_vector(@atom);
      #print "@vector\n";
    }
    else {
      my %sample_vector;
      if ($count == 1) {
        my $i;
        %sample_vector = map {$_ => ++$i} canon_vector(@atom);
      }
      my $i;
      my %vector = map {++$i => $_} canon_vector(@atom);
      @vector = @vector{keys %sample_vector}
    }
    $max = @vector;
  }
  
  if ($crop) {
    foreach (@vector) {
      print " $atom[$_]\n";
    }
  }
  else {
    foreach ($max+1 .. $N) {
      $vector[$_-1] = $_;
    }
    #print "@vector\n";
    for ($i=1;$i<=$N;$i++) {
      $ii = $vector[$i-1];
      print "$atom[$ii]";
      print "    $ii" if $prN;
      print "\n";
    }
  }
  
}

sub canon_vector {
  my @atom = @_;
  my $mol = Chemistry::Mol->new();
  for (my $i=1; $i<@atom; $i++) {
    my ($at,$x,$y,$z) = split ' ', $atom[$i];
    $mol->new_atom(symbol => $at, coords => [$x,$y,$z]);
  }
  canonicalize($mol, 'sort'=>0);
  my @vector;
  for (my $i=1; $i<@atom; $i++) {
    $vector[$i-1] = $mol->atoms($i)->attr("canon/class")
  }
  return @vector;
}

sub read_molden {
 ############################################################################
 ## Молекула -- ссылка на массив, 0-й элемент -- свойства, следующие - атомы.
 ## Свойства -- ссылка на хэш с ключами Energy, Symmetry
 ## Атом -- ссылка на массив [atom, x, y, z, ppm] (ppm может не быть).
 ##
 ## Читает xyz. Параметры - имена xyz-файлов. Если параметров нет, то <>.
 ## Возвращает массив найденных молекул.
 ############################################################################
  local @ARGV = @_ ? @_ : @ARGV;
  my $fnum = qr/-?\d+(?:\.\d+)?/;
  my $num = qr/-?\d+\.\d+(?:[eE][+-]?\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]{Title} = $line;
      ($mol[0]{Energy}) = $line =~ /(?:\s|^|=)($num)(?:\s|$)/;
      ($mol[0]{Symmetry}) = $line =~ /symm\S*\s+(\S+)/i;
      ($mol[0]{Charge}) = $line =~ /Charge\s+(-?\d+)/;
      ($mol[0]{Mult}) = $line =~ /Mult\s+(\d+)/;
      ($mol[0]{HoF}) = $line =~ /HoF\s+($fnum)/o;
      ($mol[0]{Edisp}) = $line =~ /Edisp\s+($fnum)/o;
      ($mol[0]{ZPE}) = $line =~ /ZPE\s+($fnum)/o;
      ($mol[0]{Dipole}) = $line =~ /Dipole\s+($fnum)/o;
      while ($line =~ /G\(($fnum)\)\s+($fnum)/g) {
        push @{$mol[0]{G}}, [$1,$2];
      }
      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;
  
  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 $full_title = 1;
    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/) {
        #print "  $f $mol->[0]{$f}";
        print "  $f $mol->[0]{$f}" if defined $mol->[0]{$f};
      } 
      if ($mol->[0]{G}) {
        foreach my $g (@{$mol->[0]{G}}) {
          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 $mol->[$i][4];
      print "\n";
    }
  }
  
  if ($oldfh) {
    close F;
    select $oldfh;
  }
}

sub obabel_canonical {
  # obabel неправильно работает, если в multi-xyz строка с числом атомов 
  # не начинается с пробела. Код ниже исправляет это.
  foreach my $file (@_) {
    my @mols = read_molden($file);
    write_molden(@mols, $file);
  }
  
  my @args = ($obabel_x, @_, '-O', 'obabel_canonical.xyz', '--canonical');
  warn "@args\n";
  system(@args) == 0 or die "System\n@args\nfailed: $?\n";
  
  my $m0 = (read_molden($_[0]))[0];
  my @mols = read_molden('obabel_canonical.xyz');
  unlink 'obabel_canonical.xyz';
  if ($auto == 0) {
    write_molden(@mols);
    return;
  }
  my $m1 = $mols[0];
  die "Not isomers @_\n" if @$m0 != @$m1;
  
  my @vector;
  for (my $i=1; $i<@$m0; $i++) {
    for (my $j=1; $j<@$m1; $j++) {
      next if $m1->[$j][0] ne $m0->[$i][0];
      my $sum = abs($m1->[$j][1]-$m0->[$i][1])+
                abs($m1->[$j][2]-$m0->[$i][2])+
                abs($m1->[$j][3]-$m0->[$i][3]);
      if ($sum < 1e-3) {
        push @vector, $j;
        last;
      }
    }
  }
  die "Bad vector @vector\n" if @vector != @$m0-1;
  #warn "@vector\n";
  
  foreach my $mol (@mols) {
    my @m = ($mol->[0]);
    foreach my $i (@vector) {
      push @m, $mol->[$i];
    }
    write_molden(\@m);
  }
}