#!/usr/bin/perl -ws
use strict;

#############################EDIT HERE##########################################

# Priroda
#my $priroda = './p9';
my $priroda = '/usr/local/lib/priroda/p6_32';

# Базис, если он задан без пути, будет сначала искаться 
# в текущей директории, затем в $basis_dir
my $basis_dir = '/usr/local/lib/priroda/basis';

# Директория для временных файлов Природы (на разделе, где много места)
# если задана в in-файле, то будет она
# если там не задана, то $path_tmp
# если и $path_tmp не задана, то Priroda default (/tmp)
my $path_tmp = "/disk2/tmp/$ENV{USER}";
#my $path_tmp = '.'; # Текущая директория

# Программа для обработке output-ов (с путем, если нужно).
# Если обработка output-ов этим скриптом не нужна, закоммент. след. строку
my $pri2mol = '/usr/local/bin/pri2mol';

# Нужно ли удалять input-файлы для отдельных задач
my $delete_tmp_input = 0; # Не нужно
#my $delete_tmp_input = 1; # Нужно

#############################END OF EDIT########################################

# Обозначения (суффиксы) для имен файлов в зависимости от задачи
# Можно редактировать на свой вкус
my %task_suff = (
 # Задача     Суффикс 
  optimize => 'O', 
  hessian  => 'H', 
  saddle   => 'TS', 
  nmr      => 'NMR',
  scan     => 'SCAN',
  irc0     => 'IRC0',
  irc1     => 'IRC1',
  energy   => 'E',
  gradient => 'G',
  dipole   => 'D',
  Optimize => 'Opt', 
);
# Еще в качестве задач предусмотрены IRC и TS (см. help)

our($h,$help,$np,$log);
$np ||= 1;

if ($h || $help) {
  (my $program = $0) =~ s/^.*[\/\\]//;
  print qq{
Usage: $program [-np=N] input.in [input1.in ...]
Зависимости: perl, Priroda

Прежде всего отредактируйте значения переменных в начале скрипта.

Опции: 
  -h     help
  -help  подробный help
  -np=N  Запускать на N ядрах (по умолчанию на одном)
  -log   Печать будет не на консоль, а в log-файл
  
Скрипт для последовательных запусков квантово-химической программы Природа.
Последовательность запусков определяется строкой после task= в группе \$control
input-файла. Строка состоит из последовательности задач, разделенных знаком +, 
и не должна содержать пробелов. После каждой задачи могут быть в круглых скобках
уточнения вида (метод/фунционал/базис).

Примеры (b1 и L22 -- "легкий" и "тяжелый" базисы):
task=optimize+hessian+nmr(L22)
task=hessian(b1)+optimize(b1)+optimize
task=hessian(b1)+saddle+hessian+irc
  };
  print "\n$program -help покажет больше информации\n" unless $help;
  if ($help) {
    print qq{
Не нужно указывать названия output-файлов, input-файлы можно задавать по маске.
На основании названия input-файла определяется имя и расширение output-файла.
Если input-файл кончается на '.in' или '.inp', то расширения у output-файла не
будет. В противном случае именем будет полное название input-файла, 
расширением -- '.out'. Между именем и расширением будет суффикс, свой для каждой
задачи в соответствии со следующей таблицей};
    my @arr = map { sprintf "%10s => %-4s",$_,$task_suff{$_} } sort keys %task_suff;
    for (my $i=0; $i<@arr; $i++) {
      print "\n" if $i%4 eq 0;
      print $arr[$i];
    }
    print qq{
Суффиксы в этой таблице можно отредактировать в тексте скрипта.
Суффикс, совпадающий с одним из предыдущих, будет инкрементироваться.

В основном задачи совпадают с таковыми для самой "Природы". Исключения:
irc -- это irc0+irc1 (task=irc, back=0 для irc0, back=1 для irc1)
TS -- это синоним saddle (task=optimize, saddle=1)
Optimize -- оптимизация с рестартами, см. ниже
В отличие от "Природы" регистр символов важен.

Если задача указана без уточнения (в круглых скобках), то метод/фунционал/базис
берутся из input-файла (или умолчаемые "Природы"). Уточнения м.б. вида
(метод/фунционал/базис) либо (фунционал/базис) либо (базис). Если методу не
нужен фунционал, то (метод//базис). Базис -- это файл базиса. Если он с полным
путем, то д.б. в двойных кавычках. Если не нужны разные базисные наборы для
разных атомов, удобно использовать файлы с короткими именами, по названию
базисного набора, помещенные в специальную директорию (см. начало скрипта).
См. о базисах http://limor1.nioch.nsc.ru/priroda.html

Некоторые особенности:
--Всегда используется последняя оптимизированная геометрия.
--Перед irc и saddle д.б. посчитан hessian и в нем д.б. мнимая частота 
  (или в input-файле д.б. группа с \$Energy).
--Предварительно посчитанный hessian использется и для optimize за исключением
  случая, когда он уже не актуален, типа hessian(b1)+optimize(b1)+optimize (для
  второго optimize).
--Планируется сделать особую обработку некоторых полезных сочетаний:
  scan+saddle, scan+nmr, irc+nmr.

На консоль скрипт печатает название input-файла и строку задач для него и
для каждой задачи метод/фунционал/базис, команду запуска расчета "Природой" и
краткую информацию о памяти/времени счета.

Особый случай оптимизации геометрии -- Optimize (optimize с заглавной буквы).
  Делает следующее.
  --Включает запись файла векторов (если save= в группе \$control не было 
    задано в явном виде, этот файл потом удаляется).
  --Если за заданное количество циклов (iterations= в группе \$scf) градиент SCF
    остается слишком большим ( SCF is far from  convergence) и Природа вылетает,
    то она перезапускается с текущей геометрией и последними векторами (read=1 в
    группе \$guess). Говорят, в сложных случаях помогает умеренное кол-во циклов
    (iterations=10) с последующим перезапуском.
  --Если оптимизация геометрии не сходится за заданное количество шагов 
    (steps= в группе \$optimize, 50 by default), то перестартовывает Природу 
    с последней геометрией. Это полезно, когда оптимизация не сходится даже 
    за 300 шагов, "зацикливаясь" около некоторой геометрии.
  Все output'ы Природы после перезапусков разделяются строкой из решеток и 
  записываются в один файл с суффиксом $task_suff{Optimize}. На консоль скрипт печатает 
  циклы SCF, градиенты оптимизации геометрии и некоторую информацию о своих 
  действиях. Также на консоль попадают сообщения Природы в stderr, если их не 
  перенаправить.
    };
  };
  exit;
}

die "No Priroda executable $priroda!\n" unless -x $priroda;

die "Usage: $0 [-np=N] input.in [input1.in ...]\n$0 -h for help\n" unless @ARGV;

if ($pri2mol && ! -x $pri2mol) {
  warn "Warning: No $pri2mol executable\n";
  undef $pri2mol;
}

@ARGV = glob("@ARGV");
while (@ARGV) {
  # Над каждым input-файлом будем проводить серию расчетов
  my $input = shift;
  # На основании имени input-файла определяем имя и расширение
  # Они будут постоянными, между ними будет суффикс, 
  # свой для каждого расчета в соответствии с %suff
  my ($name,$ext) = $input=~/(.*)\.inp?$/ ? ($1,'') : ($input,'.out');
  
  # Зачитываем input-файл
  my $in = read_in($input);
  
  # Печатать в log-файл (если опция -log)
  if ($log) {
    #(my $log_name = $name) =~ s/\.$in->{control}{basis}$//;
    my $log_name = $name;
    open LOG, '>>', "$log_name.log" or die "Can't write to $log_name.log: $!\n";
    select LOG;
    $| = 1;
  }
  
  print "Input fife $input\n";
  
  # Определяем директорию для временных файлов
  $in->{system}{path} ||= ($path_tmp || '/tmp');
  unless (-d $path_tmp) {
    mkdir $path_tmp or die "Can't create dir for tmp files $path_tmp: $!\n";
  }
  if ($in->{system}{path} ne '/tmp') {
    unlink grep {m/tmp\d+\.\d+/} glob("$in->{system}{path}/tmp*");
  }

  # Метод/функционал/базис
  # Если метод не предполагает функционала, то на его месте пустая строка
  
  my ($basis, $basis_file);
  if ($in->{control}{basis}) {
    $basis = put_basis($in, $in->{control}{basis});
    $basis_file = $in->{control}{basis};
  }
  
  my $method = $in->{control}{theory} || 'DFT';
  
  my $functional = $in->{dft}{functional} || 'PBE';
  
  $in->{control}{task} ||= 'gradient'; # Priroda default task
  print "Task:  $in->{control}{task}  $method/$functional/$basis\n";
  
  # Модифицируем строку task= если нужно
  my $new_task = $in->{control}{task};
  # irc заменяем на irc0+irc1
  $new_task =~ s/\birc(?=(?:[(+]|$))([^+]*)/irc0$1+irc1$1/ig;
  $new_task =~ s/\bTS(?=(?:[(+]|$))/saddle/ig;
#  # Для некоторых задач необходим предварительно посчитанный гессиан, поэтому
#  # saddle|irc заменяем на hessian+saddle|irc если hessian не был задан явно
#  if ($new_task !~ /\bhessian\b.*?\+(?:saddle|irc)/) {
#    $new_task =~ s/\b(saddle|irc[01])\b([^+]*)/hessian$2+$1$2/;
#  }
  if ($new_task ne $in->{control}{task}) {
    $in->{control}{task} = $new_task;
    print "tasks: $in->{control}{task}\n";
  }
  print "\n";
  
  # Парсим строку (д.б. без пробелов) после task=
  # Получаем массив хэшей с ключами t, b, f, m
  # Базис с путем д.б. в "
  my @task;
  foreach (split /\+/, $in->{control}{task}) {
    my ($t,$b,$f,$m);
    {
      s/^(\w+)// or last;
      $t = $1;
      s/^\(// && s/\)$// or last;
      s/("[^"]+"|[^\/]+)$// if $_;
      ($b = $1) =~ s/"//g;
      s/\/$//;
      ($f,$m) = reverse split /\//;
    }
    push @task, {t=>$t,b=>$b,f=>$f,m=>$m};
  }
  
  my %suffs;
  my ($last_geometry,$last_hessian);
  for (my $i=0; $i<@task; $i++) {
    
    $in->{control}{task} = $task[$i]{t};

    my $new_basis;
    if ($task[$i]{b}) {
      $new_basis = put_basis($in, $task[$i]{b});
    } else {
      $in->{control}{basis} = $basis_file;
      $new_basis = $basis;
    }
    die "No basis specified\n" unless $in->{control}{basis};

    $in->{control}{theory} = $task[$i]{m} || $method;
    
    $in->{dft}{functional} = $task[$i]{f} || $functional;
    delete $in->{dft}{functional} if $in->{control}{theory} !~ /dft/i;

    my $suff = $task_suff{$task[$i]{t}} || uc($task[$i]{t});

    if ($task[$i]{t} =~ /^optimize$/i && !$in->{optimize}{saddle}) {
      $suff = $task_suff{optimize};
      $in->{Energy}{eng} = $last_hessian if $last_hessian;
    }
    elsif ($task[$i]{t} eq 'saddle' or 
           $task[$i]{t} =~ /^optimize$/i && $in->{optimize}{saddle}==1) 
    {
      $in->{control}{task} = 'optimize';
      $in->{optimize}{saddle} = 1;
      $suff = $task_suff{saddle};
      $in->{Energy}{eng} = $last_hessian if $last_hessian;
    } 
    elsif ($task[$i]{t} =~ /^irc([01])$/) {
      $in->{control}{task} = 'irc';
      $in->{optimize}{back} = $1 eq 0 ? 0 : 1;
      $in->{optimize}{points} ||= 100;
      $in->{Energy}{eng} = $last_hessian if $last_hessian;
    }
    elsif ($task[$i]{t} eq 'hessian') {
      delete $in->{Energy};
      if (! $in->{control}{print} || $in->{control}{print} !~ /molden/) {
        $in->{control}{print} .= '+molden';
      }
    }
    elsif ($task[$i]{t} eq 'nmr') {
      delete $in->{Energy};
    }
    elsif ($task[$i]{t} eq 'energy') {
      delete $in->{Energy};
      #if (! $in->{control}{print} || $in->{control}{print} !~ /vectors/) {
      #  $in->{control}{print} .= '+vectors';
      #  $in->{control}{print} .= '+molden' if $in->{control}{print} !~ /molden/;
      #}
    }
    elsif ($task[$i]{t} eq 'gradient') {
      delete $in->{Energy};
    }

    if ($task[$i]{b}) {
      $suff .= "_$new_basis";
      if ($task[$i]{f}) {
        $suff .= "_".($task[$i]{f}||'');
      }
      if ($task[$i]{m}) {
        $suff .= "_$task[$i]{m}";
      }
    }
    print "task=$task[$i]{t} ",
           "$in->{control}{theory}/",
           $in->{dft}{functional}||'', '/',
           "$new_basis\n";

    while (exists $suffs{$suff}) {
      $suff =~ s/([A-Z]+)(\d*)/$1.(($2||0)+1)/e;
    }
    $suffs{$suff} = 1;
    $suff = ".$suff";

    if ($i > 0) {
      $in->{molecule}{MOL} = $last_geometry if $last_geometry;
      #$in->{Energy}{eng} = $last_hessian if $last_hessian;
      if ($suff =~ /$task_suff{saddle}|$task_suff{irc0}|$task_suff{irc1}/) {
        die "No hessian\n" unless $in->{Energy}{eng};
        die "No imaginary freq.\n" if $in->{Energy}{eng} !~ /Imaginary/i;
      }
      my $prev_tasks = join '+', map {$_->{t}} @task[0..$i-1];
      delete $in->{Energy} if 
        $task[$i]{t} eq 'hessian' ||
        $task[$i]{t} =~ /^(?:optimize|saddle)$/i &&
        $prev_tasks =~ /hessian.+(?:optimize|saddle)/ &&
        $' !~ /hessian/;
    }

    my $tmp_input = "$name$suff$ext.in";
    my $output = $name.$suff.$ext;

    #print "Run Priroda:\n";
    if ($task[$i]{t} eq 'Optimize') {
      smart_optimize($in, $tmp_input, $output);
    }
    else {
      print_in($in, $tmp_input);
      my @args = ($priroda, '-np', $np, $tmp_input, $output);
      print "`@args`";
      system(@args) == 0 or die "\nSystem\n@args\nfailed: $?\n";
      print "\n";
    }
    
    my ($MOL,$eng,$info) = parse_output($output);
    print "$info\n";
    if    ($in->{control}{task} eq 'optimize') {
      die "Optimized geometry not found in $output\n" unless $MOL;
      $last_geometry = $MOL;
      delete $in->{optimize}{saddle};
    }
    elsif ($in->{control}{task} eq 'hessian') {
      die "Hessian not found in $output\n" unless $eng;
      $last_hessian = $eng;
    }
    elsif ($in->{control}{task} eq 'irc') {
      delete $in->{optimize}{back};
    }
    
    unlink $tmp_input if $delete_tmp_input;
    system($pri2mol, $output) if $pri2mol;
  }
  print(('-' x 60), "\n") if @ARGV;
  
  close LOG;
}

sub smart_optimize {
  my ($in, $tmp_input, $output) = @_;
  
  my $delete_vectors;
  unless ($in->{control}{save}) {
    $delete_vectors = 1;
    $in->{control}{save} = "$in->{system}{path}/$output.vec";
  }
  $in->{guess}{read} = 0;

  print_in($in, $tmp_input);

  open OUT, '>>', $output or die "Can't open $output: $!\n";
  select((select(OUT), $| = 1)[0]);
  print OUT "Priroda run under $0\n";
  my $count;
  my ($count_scf, $count_gradients);
  while (1) {
    my ($no_convergence,$curr_mol,$OPTIMIZATION_CONVERGED, $error);
  #  my ($no_convergence,$new_optimization,$curr_mol,$OPTIMIZATION_CONVERGED);
    #print "\n----------- Run $count_scf ------------\n" if $count;
    print OUT '#'x 80, "\n";
    my @args = ($priroda, '-np', $np, $tmp_input);
    #my $pid = open(PRI, '-|', @args) or die "Can't run @args: $!\n";
    #print "@args      PID $pid\n";
    open(PRI, '-|', @args) or die "Can't run @args: $!\n";
    #print "@args\n";
    select((select(PRI), $| = 1)[0]);
    while (<PRI>) {
      print OUT;
      print if / it.       energy       energy change    gradient     time/;
      if (m/^\s+\d+(?:\s*-?\d+\.\d+){4}\s*$/) {
        $count_scf++;
        print;
        #print (join ' ', stat($in->{control}{save}), "\n");
      }
      
      print if m/^ converged!/;
      print if m/^ no convergence!/;
      if (m/^ SCF is far from convergence/) {
        #print;
        $no_convergence = 1;
      }
      
      if (m/error/i) {
        print;
        $error = 1;
      }

      if (m/^ max\. gradient/){
        print;
        $count_gradients++;
      }

      if (m/^mol>\s*\$molecule/..m/^mol>\s*\$end/) {
        undef $curr_mol if m/^mol>\s*\$molecule/;
        s/^mol>//;
        $curr_mol .= $_;
      }

      if (m/^ OPTIMIZATION CONVERGED/) {
        $OPTIMIZATION_CONVERGED = 1;
      }
    }
    close PRI;

    if ($error) {
      last;
    }

    if ($curr_mol && $curr_mol =~ /(cartesian|z-matrix)(.+)\$end/is) {
      #print $2;
      $in->{molecule}{$1} = $2;
    }

    if ($OPTIMIZATION_CONVERGED) {
      print " OPTIMIZATION CONVERGED\n";
      print $count_scf||0, ' SCFs, ', $count_gradients||0, " gradients\n";
      unlink $tmp_input;
      last;
    }

    if ($no_convergence) {
      print " SCF is far from convergence. Start with current geometry and last vectors\n\n";
      $count++;
      if ($count == 1) {
        $in->{guess}{read} = 1;
        print_in($in, $tmp_input);
      }
      sleep 1;
      next;
    }

    print " OPTIMIZATION NOT CONVERGED. Start with last geometry\n";
    $count = 0;
    $in->{guess}{read} = 0;
    delete $in->{Energy} if exists $in->{Energy};
    print_in($in, $tmp_input);
  }

  close OUT;

  unlink $in->{control}{save} if $delete_vectors;

}

sub read_in {
# Читает input Природы. Параметр - file.in.
# Возвращает ссылку на хэш.
# Устройство этого хэше вполне очевидно: 
#   $system mem=128 $end      $in->{system}{memory} = 128;
# за исключением слов со скобками:
#   $guess na(85)=86,85 $end  $in->{guess}{na} = ['(85)', '86,85'];
# Группы $molecule...$end и $Energy...$end - особые случае (см. код)
  
  # Вспомогалельные шаблоны
  my $d = qr/\d+/;         # Целое число
  my $s = qr/"[^"]*"|\S+/; # Строка
  my $f = qr/(?:[+-]?)(?=\d|\.\d)\d*(?:\.\d*)?(?:[Ee]([+-]?\d+))?/; # Число C
  # Хэш ключевых слов (по группам) и шаблоны значений
  my %IN = (
   system => [
    qr/( memory      ) (\(.*?\))? \s* = \s* ( $d             ) /xoi,
    qr/( disk        ) (\(.*?\))? \s* = \s* ( $d             ) /xoi,
    qr/( path        ) (\(.*?\))? \s* = \s* ( $s             ) /xoi,
   ],
   control => [
    qr/( task        ) (\(.*?\))? \s* = \s* ( $s             ) /xoi,
    qr/( theory      ) (\(.*?\))? \s* = \s* ( $s             ) /xoi,
    qr/( state       ) (\(.*?\))? \s* = \s* ( $d(?:,$d)?     ) /xoi,
    qr/( basis       ) (\(.*?\))? \s* = \s* ( $s             ) /xoi,
    qr/( four        ) (\(.*?\))? \s* = \s* ( $d             ) /xoi,
    qr/( nucleus     ) (\(.*?\))? \s* = \s* ( $s             ) /xoi,
    qr/( print       ) (\(.*?\))? \s* = \s* ( $s             ) /xoi,
    qr/( save        ) (\(.*?\))? \s* = \s* ( $s             ) /xoi,
   ],
   guess => [
    qr/( read        ) (\(.*?\))? \s* = \s* ( [01]           ) /xoi,
    qr/( na          ) (\(.*?\))? \s* = \s* ( $d(?:[,\s]+$d)*) /xoi,
    qr/( nb          ) (\(.*?\))? \s* = \s* ( $d(?:[,\s]+$d)*) /xoi,
   ],
   dft => [
    qr/( functional  ) (\(.*?\))? \s* = \s* ( $s             ) /xoi,
   ],
   scf => [
    qr/( restrict    ) (\(.*?\))? \s* = \s* ( [01]           ) /xoi,
    qr/( convergence ) (\(.*?\))? \s* = \s* ( $f(?:,$f)?     ) /xoi,
    qr/( iterations  ) (\(.*?\))? \s* = \s* ( $d(?:,$d)?     ) /xoi,
    qr/( procedure   ) (\(.*?\))? \s* = \s* ( $s             ) /xoi,
    qr/( d1small     ) (\(.*?\))? \s* = \s* ( [01]           ) /xoi,
    qr/( core        ) (\(.*?\))? \s* = \s* ( $d             ) /xoi,
   ],
   atoms => [
    qr/( core        ) (\(.*?\))? \s* = \s* ( $d(?:[,\s]+$d)*) /xoi,
    qr/( mcore       ) (\(.*?\))? \s* = \s* ( $d             ) /xoi,
   ],
   grid => [
    qr/( accuracy    ) (\(.*?\))? \s* = \s* ( $f             ) /xoi,
   ],
   optimize => [
    qr/( saddle      ) (\(.*?\))? \s* = \s* ( [01]           ) /xoi,
    qr/( tolerance   ) (\(.*?\))? \s* = \s* ( $f(?:,$f)?     ) /xoi,
    qr/( trust       ) (\(.*?\))? \s* = \s* ( $f(?:,$f){0,2} ) /xoi,
    qr/( follow      ) (\(.*?\))? \s* = \s* ( $d             ) /xoi,
    qr/( steps       ) (\(.*?\))? \s* = \s* ( $d             ) /xoi,
    qr/( cartesian   ) (\(.*?\))? \s* = \s* ( [01]           ) /xoi,
    qr/( back        ) (\(.*?\))? \s* = \s* ( [01]           ) /xoi,
    qr/( points      ) (\(.*?\))? \s* = \s* ( $d             ) /xoi,
    qr/( radius      ) (\(.*?\))? \s* = \s* ( $f             ) /xoi,
    qr/( coordinated ) (\(.*?\))? \s* = \s* ( $d(?:[,\s]+$d)*) /xoi,
    qr/( fix         ) (\(.*?\))? \s* = \s* ( $d(?:[,\s]+$d)*) /xoi,
    qr/( value       ) (\(.*?\))? \s* = \s* ( $f(?:[,\s]+$f)*) /xoi,
    qr/( gradient    ) (\(.*?\))? \s* = \s* ( $d             ) /xoi,
    qr/( print       ) (\(.*?\))? \s* = \s* ( [01]           ) /xoi,
   ],
   thermo => [
    qr/( sigma       ) (\(.*?\))? \s* = \s* ( $d)              /xoi,
    qr/( temperature ) (\(.*?\))? \s* = \s* ( $f(?:[,\s]+$f)*) /xoi,
   ],
   d2edr2 => [
    qr/( length      ) (\(.*?\))? \s* = \s* ( $f             ) /xoi,
    qr/( displace    ) (\(.*?\))? \s* = \s* ( $d             ) /xoi,
    qr/( steps       ) (\(.*?\))? \s* = \s* ( $d(?:,$d)?     ) /xoi,
   ],
   nmr => [
    qr/( standard    ) (\(.*?\))? \s* = \s* ( $f(?:[,\s]+$f)*) /xoi,
    qr/( points      ) (\(.*?\))? \s* = \s* (         $f[,\s]+$f[,\s]+$f
                                 (?:[,\s]$f[,\s]+$f[,\s]+$f)* )/xoi,
   ],
   molecule => [
    qr/( charge      ) (\(.*?\))? \s* = \s* ( [+-]?$d        ) /xoi, 
    qr/( mult        ) (\(.*?\))? \s* = \s* ( $d             ) /xoi, 
    qr/(cartesian|z-matrix)()(.+)/is,
   ],
  
   Energy => [
    qr/( E           ) (\(.*?\))? \s* =     ( [-+\d.e\s]+     ) /xo, 
    qr/( D           ) (\(.*?\))? \s* =     ( [-+\d.e\s]+     ) /xo, 
    qr/( A           ) (\(.*?\))? \s* =     ( [-+\d.e\s]+     ) /xo, 
    qr/( G           ) (\(.*?\))? \s* =     ( [-+\d.e\s]+     ) /xo, 
    qr/( Dx          ) (\(.*?\))? \s* =     ( [-+\d.e\s]+     ) /xo, 
    qr/( H           ) (\(.*?\))? \s* =     ( [-+\d.e\s]+     ) /xo, 
   ],
  
  );
  
  # Ссылка на вспомогательную функцию.
  # Принимает сокращенное слово (в начале) и полные слова плоским списком
  # Возвращает полное слово, соответствующее сокращенному
  my $extended_word = sub {
    my $x = shift;
    my $min_rest = 1e10;
    my $a;
    foreach (@_) {
      m/^$x(.*)/ or next;
      my $l = length $1;
      if ($l < $min_rest) {
        $min_rest = $l;
        $a = $_;
      }
    }
    return defined $a ? $a : $x;
  };
  
  # Открываем файл и читаем его в одну строку
  open F, '<', $_[0] or die "Can't open $_[0] : $!\n";
  local $/;
  my $in = <F>;
  close F;
  
  my %in; # Хэш ключевых слов из файла и их значений
  # Разбиваем строку по группам
  while ($in =~ /\$(\w+)(.*?)\$end/isg) {
    my ($k,$v) = ($1,$2);
    $in{$k}{rest} = $v; # В rest - исходная строка для группы
    warn("Unknown group \$$k ... \$end\n"),next if ! exists $IN{$k};
    
    # Из хэша %IN извлекаем полные ключевые слова для данной группы
    my @words = map {m/\(.*?\(\s*(\w+)/} @{$IN{$k}};
    # и заменяем сокращенные слова на полные
    $v =~ s/(\w+)(?=(?:\(.*?\))?\s*=)/$extended_word->($1,@words)/eg;
    $in{$k}{rest} = $v; # В rest - строка для группы с полными словами
    
    # По всем шаблонам для данной группы
    foreach my $pat ( @{$IN{$k}} ) {
      # сопоставляем строку с шаблоном
      $v =~ $pat or next;
      # $1 - слово; $2 - параметры слова (в скобках); $3 - значение
      $in{$k}{$1} = $2 ? [$2,$3] : $3;
      # удаляем это все из rest
      $in{$k}{rest} =~ s/$pat//;
    }
    $in{$k}{rest} =~ s/^\s+//;
    $in{$k}{rest} =~ s/\s+$//;
    #delete $in{$k}{rest} if $in{$k}{rest} eq '';
  }
  $in{molecule}{charge} = 0 if ! exists $in{molecule}{charge};
  $in{molecule}{mult} = 1 if ! exists $in{molecule}{mult};
  return \%in;
}

sub print_in {
  my @order = (
   [qw( system   memory disk path                                       )],
   [qw( control  task theory state basis four nucleus print save        )],
   [qw( guess    read na nb                                             )],
   [qw( dft      functional                                             )],
   [qw( scf      restrict convergence iterations procedure d1small core )],
   [qw( atoms    core mcore                                             )],
   [qw( grid     accuracy                                               )],
   [qw( optimize saddle tolerance trust follow steps cartesian back  
                 points radius coordinated fix value gradient print     )],
   [qw( thermo   sigma temperature                                      )],
   [qw( d2edr2   length displace steps                                  )],
   [qw( nmr      standard points                                        )],
   #[qw( molecule charge mult cartesian z-matrix                         )],
   #[qw( Energy   E D A G Dx H                                           )],
  );
  my %in = %{$_[0]};
  my $stdout;
  if (defined $_[1]) {
    open F, '>', $_[1] or die "Can't open $_[1]: $!\n";
    $stdout = select F;
  }
  
  foreach (@order) {
    my @words = @$_;
    my $group = shift @words;
    next if ! exists $in{$group};
    my @a;
    push @a, "\$$group";
    foreach (@words) {
      next if ! exists $in{$group}{$_};
      push @a, $_. (ref($in{$group}{$_}) ? 
                   "$in{$group}{$_}[0]=$in{$group}{$_}[1]" : 
                   "=$in{$group}{$_}");
    }
    push @a, "$in{$group}{rest}" if $in{$group}{rest};
    push @a, "\$end";
    delete $in{$group};
    my $str = "@a";
    next if $str =~ /^\$$group\s*\$end$/;
    if (length $str < 60) {
      print "$str\n";
    } else {
      $_=" $_" for @a[1..$#a-1];
      print join("\n", @a), "\n";
    }
  }
  foreach (grep {!/molecule|Energy/} sort keys %in) {
    print "\$$_$in{$_}{rest}\$end\n" if $in{$_}{rest};
  }
  
  if ($in{molecule}) {
    if (exists $in{molecule}{MOL}) {
      print $in{molecule}{MOL};
    } else {
      print "\$molecule";
      print " charge=$in{molecule}{charge}" if $in{molecule}{charge} != 0;
      print " mult=$in{molecule}{mult}" if $in{molecule}{mult} != 1;
      print "\n";
      foreach (qw(cartesian z-matrix)) {
        next if ! exists $in{molecule}{$_};
        print " $_$in{molecule}{$_}";
      }
      print "\$end\n";
    }
  }
  
  if ($in{Energy}) {
    if (exists $in{Energy}{eng}) {
      #print "exists $in{Energy}{eng}\n";
      
      print $in{Energy}{eng};
    } else {
      print "\$Energy\n";
      foreach (qw(E D A G Dx H)) {
        next if ! exists $in{Energy}{$_};
        $in{Energy}{$_} =~ s/\s+$//;
        print " $_=$in{Energy}{$_}\n";
      }
      print "\$end\n";
    }
  }
  
  if (defined $_[1]) {
    select $stdout;
    close F;
  }
}

sub get_optimized_geometry {
  my $file = shift;
  open F, '<', $file or die "Can't open $file: $!\n";
  my $MOL;
  while (<F>) {
    undef $MOL if m/^MOL>/;
    while (m/^MOL>/) {
      $MOL .= $';
      last if eof(F);
      $_ = <F>;
    }
  }
  close F;
  return $MOL;
}

sub get_hessian {
  my $file = shift;
  open F, '<', $file or die "Can't open $file: $!\n";
  my ($eng,$im_freq);
  while (<F>) {
    undef $eng if m/^eng>/;
    while (m/^eng>/) {
      $eng .= $';
      last if eof(F);
      $_ = <F>;
    }
    if (!$im_freq && m/^\s*(freq\.(?:\s+\d+\.\d+\s+i)+)/) {
      $im_freq = $1;
    }
  }
  close F;
  return undef unless $eng;
  $eng = " Imaginary $im_freq\n$eng" if $im_freq;
  return $eng=~/H=/ ? $eng : undef;
}

sub parse_output {
  my $file = shift;
  open F, '<', $file or die "Can't open $file: $!\n";
  my ($MOL,$eng,$info); # For return
  my ($im_freq,$mem,$disk,$time) = ('',0,0,0);
  while (<F>) {
    undef $MOL if m/^MOL>/;
    while (m/^MOL>/) {
      $MOL .= $';
      last if eof(F);
      $_ = <F>;
    }
    undef $eng if m/^eng>/;
    while (m/^eng>/) {
      $eng .= $';
      last if eof(F);
      $_ = <F>;
    }
    if (m/^ ===========/) {
      while (<F>) {
        $mem += $1 if m/Memory used =\s*(\d+)/;
        $disk += $1 if m/Disk  used =\s*(\d+)/;
        $time += $1 if m/REAL time =\s*(\d+)/;
        last if m/^ ratio/;
      }
    }
    if (!$im_freq && m/^\s*(freq\.(?:\s+\d+\.\d+\s+i)+)/) {
      $im_freq = $1;
    }
  }
  close F;
  if ($eng && $eng=~/H=/) {
    $eng = $im_freq ? " Imaginary $im_freq\n$eng" : " No imaginary freq.\n$eng";
  }
  $info .= sprintf(" memory %.3f M  ", $mem/1024)       if $mem;
  $info .= sprintf(" disk %.3f G  ",   $disk/1024/1024) if $disk;
  $info .= sprintf(" time %.2f h  ",   $time/3600)      if $time;
  $info = $info ? "$info\n" : '';
  return($MOL, $eng, $info);
}

sub put_basis {
  my ($in, $basis) = @_;

  my $b;
  ($b = $basis) =~ s/.*[\/\\]//;
  if (! $in->{control}{four} and $b =~ /^rL/ || $b eq 'basis4.in') {
    print "Set four=1\n"; 
    $in->{control}{four} = 1;
  }

  if (-e $basis) {
    $in->{control}{basis} = $basis;
  }
  elsif (-e "$basis_dir/$basis") {
    $in->{control}{basis} = "$basis_dir/$basis";
  }
  else {
    die "No basis $basis in current dir and in $basis_dir/\n";
  }
  return $b;
}